State University of New York Report 


.1 i 


| } 


t- 3 

Pi 


i M 


Large Eddy Simulations and Direct 
Numerical Simulations of High 
Speed Turbulent Reacting Flows 

by 


P. Givi, S.H. Frankel, V. Adumitroaie, 
G. Sabini, and C.K. Madnia 

yAG>l-/l2Z. 


Department of Mechanical and Aerospace Engineering 
State University of New York 
Buffalo, New York 14260-4400 


Semi-Annual Report Submitted to 
NASA Langley Research Center 

Summary of Activities Supported Under Grant NAG 1-1122 

for the Period 



August 1, 1992 - January 31, 1993 


Contents 

1 Recent Progress 2 

2 Recent Publications 5 

3 Appendix I 7 

4 Appendix II 8 

5 Appendix III 9 


Large Eddy Simulations and Direct 
Numerical Simulations of High Speed 
Turbulent Reacting Flows 


P. Givi, S.H. Frankel, V. Adumitroaie, G. Sabini, and C.K. Madnia 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, New York 14260 


Abstract 

The primary objective of this research is to extend current capabilities of Large Eddy 
Simulations (LES) and Direct Numerical Simulations (DNS) for the computational 
analyses of high speed reacting flows. Our efforts in the first two years of this research 
have been concentrated on a priori investigations of single-point Probability Density 
Function (PDF) methods for providing subgrid closures in reacting turbulent flows. In 
the efforts initiated in the third year, our primary focus has been on performing actual 
LES by means of PDF methods. The approach is based on assumed PDF methods 
and we have performed extensive analysis of turbulent reacting flows by means of LES. 
This includes simulations of both three-dimensional (3D) isotropic compressible flows 
and two-dimensional reacting planar mixing layers. In addition to these LES analyses, 
some work is in progress to assess the extent of validity of our assumed PDF methods. 
This assessment is done by making detailed companions with recent laboratory data 
in predicting the rate of reactant conversion in parallel reacting shear flows. 

This report provides a summary of our achievements for the first six months of 
the third year of this program. This research has been supported by NASA Langley 
Research Center under Grant NAG- 1-1 122. Dr. J. Philip Drummond, Theoretical Flow 
Physics Branch (TFPB), Mail Stop 156, Tel: 804-864-2298 is the Technical Monitor of 
this Grant 
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1 Recent Progress 


In the first two years of this research program, we have been speculating the feasibility of the 
approach based on probability density function (PDF) methods as a means of accounting for 
the effects of subgrid scalar fluctuations in large eddy simulations (LES) of turbulent reacting 
flows. Within this period, we have devoted substantial effort on extending the state-of-the- 
art on PDF modeling, and on making use of recent developments for the purpose of our 
investigation. Fortunately, we have been able to make noteworthy progress in understanding 
the advantages and drawbacks of all of the recently available PDF methods. Therefore, in 
this third year of our activities, we have initiated a strong program in making use of PDF 
methods for the purpose of LES. In doing so, based on our experience in the first two years, 
we have initiated a step-by-step procedure which we shall describe below. The advantage of 
this procedure is the fact that the results of our investigation are also helpful in approaches 
based on conventional Reynolds averaging methods. 

It is well-know that there are two ways by which the PDF methods can be invoked for 
modeling of scalar fluctuations in a statistical description of turbulent reacting flows [1, 2]: 
(1) Assumed PDF, (2) PDF obtained by a modeled transport equation. By fluctuations 
we mean those in either typical Reynolds averaging or within the subgrid in LES. Each of 
these two PDF methods have certain advantages and drawbacks, as indicated in our previous 
yearly reports and as also indicated in a number of previous publications, e.g. [3, 4, 5]. 

A noteworthy progress in PDF methods is the recent development of the Amplitude Mapping 
Closure (AMC) by Kraichnan and co-workers [6, 7]. In a number of standard test cases, it 
has been demonstrated that this closure is somewhat more superior in comparisons with 
previous PDF methods based on the so called Coalescence/Dispersion (C/D) models [8]. 
Based on this development, at the initial phase of this research our intention was to make 
use of AMC for the purpose of subgrid closures. For that, we initiated a research program 
for the purpose of better understanding the properties of this closure. Based on our work 
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to date, the results of which are reported in [9, 10], it is shown that the approach provides 
some advantages over previous models for some “simple” cases. However, it does also have 
certain drawbacks which limits its use for “practical” applications. In our recent paper [10] 
(included as Appendix I), a lengthy discussion is presented of all the properties of these 
methods. Here, we shall summarize some of the difficulties associated with the use of AMC 
for PDF based LES: 

1. The use of the approach is practical for closure only at the single-point level. 

2. In its present form, the model cannot take the effects of mixing on migration of scalar 
bounds into account. 

3. The computational implementation of the model for multi-scalar mixing and reaction is 
virtually impossible. 

4. For the cases which the model has been validated, other simpler PDF methods; partic- 
ularly those belonging to the Pearson Family (PF) [11] or Johnson Edgeworth Translation 
(JET) [12, 13] of distribution have been satisfactory. 

Considering points (1) and (2), it is safe to conclude that the AMC shares the same limita- 
tions as those in assumed PDF methods and in C/D closures. Point (3) advocates the use 
of other PDF methods until developments of more reliable numerical methods. Based on 
these three points, and also considering point (4) we decided to “discontinue” further basic 
research on assessing pro’s and con’s of PDF methods, in favor of extending LES technology 
by means of PDF methods. This is deemed acceptable, at least at the present stage, since: 

(i) Most current PDF methods can be only incorporated at the single-point level. Therefore, 
closure models for the second order moments ( e.g . variances and covariance) have to be 
provided by other means. 

(ii) The flexibilities offered by some of the assumed probability distributions, e.g. members 
of PF or JET is very pleasing for accounting the effects of species fluctuations in a stochastic 
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sense. 


(iii) While in Reynolds averaging of reacting flow, the modeling of the PDF is very crucial 
[14], we speculate that in LES of such flows it is more crucial to have the second order 
moments of the scalar field predicted correctly. 

Based on these points, we have decided to perform extensive LES analyses of homogeneous- 
isotropic and spatially developing reacting mixing layers. The results of our efforts to-date 
are described in detail in Appendix II. The summary of our findings is: 

(a) We feel that the approach based on a hybrid second-order closure and assumed PDF 
methods provides the most promising and practical means of performing LES of reacting 
flows, at least presently. The appropriate form of the closure that can be used is based on a 
k — A — g model. Here, k = kinetic energy of velocity fluctuations within the subgrid, A = 
the width of the filter, and g = the “covariance” vector. 

(b) It is impossible to make use of zero-order models (equivalent of Prandtl’s mixing length 
hypothesis) in LES. The implementation of this model results in wrong predictions. 

(c) The second order methodology is significantly different from that used in conventional 
Reynolds averaging procedures. The primary difference is due to modeling of convective 
transport and also the “diminished” degree of freedom because of the lack of a transport 
equation for A. 

(d) After significant analyses (see Appendix II), we now have a reasonable k equation. This 
equation works well if a combination of both the “strain rate tensor” and the “rotation 
tensor” are used in modeling the turbulent viscosity. 

(e) In simple reacting systems, assumed PDF’s of PF and JET [10] provide a reasonable 
means of accounting for the scalar fluctuations within the subgrid. For example, in the 
chemistry of A + B — ► Product, the Logit-Normal density and the Beta density of first 
kind are reasonable for the chemistry under equilibrium conditions. Under non-equilibrium 
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conditions, the use of Dirichlet distribution can be justified. 

(f) The model described above can be used effectively if k and g are predicted accurately. 
The usual practice is to use gradient diffusion modeling for the effects of convective transport 
in the g equation. While this model has been somewhat satisfactory in Reynolds averaging 
procedures, it does not work in LES (see Appendix II). In fact, even in a simple non- 
reacting flow, the scalar g (not in the vector form) cannot be precited well. This is expected 
considering the different form of k equation which had to be used. Nevertheless, if this 
parameter is modeled accurately the LES procedure can be implemented very effectively. 

We do not yet know how to modify the g equation in a form appropriate for LES. This issue 
is currently being investigated. 

The next step in our work in this direction is to replace the assumed PDF method with the 
approach based on a transport equation for the PDF. For that we plan to make use of C/D 
models for modeling of the molecular mixing term. Our purpose for doing so is due to two 
factors: (1) C/D closure can be incorporated into the numerical procedure rather easily, and 
(2) The problem associated with the lack of boundary migration is less severe in C/D in 
comparison to AMC. 

The summary of all our achievements are highlighted in the appendices. These appendices 
are self explanatory. Therefore, there is no need for further elaboration here. 


2 Recent Publications 

Within the first half of our activities in the third year, we have produced the following 
publications: A journal article to appear Comb. Set. and Tech. [10] (Appendix I), a paper 
to be presented at the LES conference at the ASME meeting in June 1993 (Appendix II), 
and an article to be submitted for publication shortly (Appendix III). 


5 


References 


[1] O’Brien, E. E., The Probability Density Function (PDF) Approach to Reacting Tur- 
bulent Flows, in Libby, P. A. and Williams, F. A., editors, Turbulent Reacting Flows , 
chapter 5, pp. 185-218, Springer- Verlag, Heidelberg, 1980. 

[2] Pope, S. B., PDF Methods for Turbulent Reacting Flows, Prog. Energy Combust. Sci., 
11:119-192 (1985). 

[3] Givi, P., Model Free Simulations of Turbulent Reactive Flows, Prog. Energy Combust. 
Sci., 15:1-107 (1989). 

[4] Kollmann, W., The pdf Approach to Turbulent Flow, Theoret. Comput. Fluid Dynam- 
ics, 1:249-285 (1990). 

[5] Pope, S. B., Computations of Turbulent Combustion: Progress and Challenges, in Pro- 
ceedings of 23rd Symp. (Int.) on Combustion, pp. 591-612, The Combustion Institute, 
Pittsburgh, PA, 1990. 

[6] Kraichnan, R. H., Closures for Probability Distributions, Bull. Amer. Phys. Soc., 
34:2298 (1989). 

[7] Chen, H., Chen, S., and Kraichnan, R. H., Probability Distribution of a Stochastically 
Advected Scalar Field, Phys. Rev. Lett., 63(24) :2657-2660 (1989). 

[8] Pope, S. B., Mapping Closures for Turbulent Mixing and Reaction, Theoret. Comput. 
Fluid Dynamics, 2:255-270 (1991). 

[9] Madnia, C. K., Frankel, S. H., and Givi, P., Reactant Conversion in Homogeneous 
Turbulence: Mathematical Modeling, Computational Validations and Practical Appli- 
cations, Theoret. Comput. Fluid Dynamics, 4:79-93 (1992). 

[10] Miller, R. S., Frankel, S. H., Madnia, C. K., and Givi, P., Johnson- Edgeworth Transla- 
tion for Probability Modeling of Binary Mixing ip Turbulent Flows, Combust. Sci. and 
Tech., (1993), In press. 

[11] Pearson, K., Contributions to the Mathematical Theory of Evolution: II. Skew Varia- 
tions in Homogeneous Material, Philos. Trans, of the Royal Soc. of London, Series A., 
186:343-414 (1895). 

[12] Johnson, N. L., Systems of Frequency Curves Generated by Methods of Translation, 
Biometrika, 36:149-176 (1949a). 

[13] Edgeworth, F. Y., On the Representation of Statistical Frequency by a Series, Journal 
of the Royal Statistical Society, Series A., 70:102-106 (1907). 

[14] Libby, P. A. and Williams, F. A., editors, Turbulent Reacting Flows, Topics in Applied 
Physics, Vol. 44, Springer- Verlag, Heidelberg, 1980. 


6 



3 Appendix I 


This appendix provides a summary of our work on assumed PDF methods. In addition to 
explaining many important mathematical characteristics of all the available PDF methods, 
this paper also provides a comparative assessment of these closure for practical applications. 
The major point of the report, as relevant to our LES analyses, is the comparison of the 
properties of PF and JET generated PDF’s with those of AMC. 

The work reported here makes a logical extension of our previous work reported in [9]. 
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Abstract — A family of Probability Density Functions (PDFs) generated by Johnson-Edgeworth Translation 
(JET) is used for statistical modeling of the mixing of an initially binary scalar in isotropic turbulence. 
The frequencies obtained by this translation are shown to satisfy some of the characteristics of the PDFs 
generated by the Amplitude Mapping Closure (AMC) (Kraichnan, 1989; Chen et al , 1989). In fact, the 
solution obtained by one of the members of this family is shown to be identical to that developed by the AMC 
(Pope, 1991). Due to this similarity and due to the demonstrated capabilities of the AMC, a justification is 
provided for the use of other members of JET frequencies for the modeling of the binary mixing problem. 
This similarity also furnishes the reasoning for the applicability of the Pearson Family (PF) of frequencies 
for modeling of the same phenomena. The mathematical requirements associated with the applications of 
JET in the modeling of the binary mixing problem are provided, and all the results are compared with 
data generated by Direct Numerical Simulations (DNS). These comparisons indicate that the Logit-Normal 
frequency portrays some subtle features of the mixing problem better than the other closures. However, none 
of the models considered (JET, AMC, and PF) are capable of predicting the evolution of the conditional 
expected dissipation and/or the conditional expected diffusion of the scalar field in accordance with DNS. It 
is demonstrated that this is due to the incapability of the models to account for the variations of the scalar 
bounds as the mixing proceeds. A remedy is suggested for overcoming this problem which can be useful in 
probability modeling of turbulent mixing, especially when accompanied by chemical reactions. While in the 
context of a single-point description the evolution of the scalar bounds cannot be predicted, the qualitative 
analytical-computational results portray a physically plausible behavior. 


1 INTRODUCTION 

The problem of binary mixing in turbulent flows has been the subject of widespread 
investigations over the past two decades (Dopazo, 1973; Pope, 1979; Pope, 1985; Pope, 
1990; Givi, 1989; Kollmann, 1990). This problem has been particularly useful in assessing 
the extent of validity of the closures developed within this period for modeling of turbulent 
mixing by Probability Density Function (PDF) methods (Dopazo, 1973; Pope, 1976; Pope, 
1979; Janicka et al 1979; Pope, 1982; Kosaly and Givi, 1987; Norris and Pope, 1991). 
Usually the problem is considered in the setting of a spatially homogeneous turbulent flow 
in which the temporal evolution of the PDF is considered. In this setting, development 
of a closure which can accurately predict the evolution of the PDF has been the main 
objective of these investigations (for recent reviews see Pope (1990); Kollmann (1990); 
Givi (1989)). 

Computational experiments based on Direct Numerical Simulations (DNS) have proven 
very useful in evaluating the performance of new closures (Givi, 1989; Pope, 1990). The 
binary mixing problem is well-suited for DNS investigation, and current computational 
capabilities allow consideration of flows at sufficiently large Reynolds numbers in which 
the behavior of the models can be assessed (Eswaran and Pope, 1988; Givi and McMurtry, 
1988; McMurtry and Givi, 1989; Madnia and Givi, 1992). The results of all the previous 
work on DNS of the binary mixing problem portray a clear picture of the PDF evolution, 
at least at the single-point level. A successful closure is one which can predict all the 
stages of mixing, as depicted by DNS, from an initially binary state (total segregated) to 
a final mixed condition. 
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Amongst the models developed in the literature, the recent Amplitude Mapping Closure 
(AMC) (Kraichnan, 1989; Chen etal, 1989; Pope, 1991) has proven effective in producing 
a physically correct PDF evolution. In the application of this model to the problem of 
binary mixing it has been demonstrated that the closure is capable of approximating a 
reasonably correct evolution at all stages of mixing (Pope, 1991). Namely, the evolution 
from an initial double delta PDF to an asymptotic Gaussian distribution. This is a trend 
which has been observed in DNS (Eswaran and Pope, 1988; Givi and McMurtiy, 1988; 
McMurtry and Givi, 1989; Madnia and Givi, 1992) and also corroborated by experimental 
investigations (Miyawaki et a/., 1974). However, it is shown by Gao (1991); O’Brien and 
Jiang (1991) that the PDF adopts an asymptotic Gaussian-like distribution only near the 
mean scalar value, and the conditional expected dissipation does not correspond to that 
of a Gaussian field everywhere within the composition domain. 

Our first objective in this work is to present another means by which the AMC can 
be viewed. It is demonstrated that in the binary mixing problem, this closure can be 
considered as a member of the family of frequencies generated by the method of Johnson- 
Edgeworth Translation (JET). In fact, it is shown that the result produced by the AMC is 
identical to that generated by one of the members of this translation. With this observation, 
a justification for employing other simpler “assumed” frequencies is provided. Our second 
objective is to make a detailed examination of the conditional expected dissipation and 
the conditional expected diffusion of the scalar variable as predicted by the closures. 
This examination provides an effective means of demonstrating the deficiencies of these 
models in reproducing the correct physical behavior as depicted by DNS results. With 
the development of analytical relations for some of these closures, a remedy is suggested 
for overcoming the model deficiencies. 

1. 1 Outline 

In the next section, the problem of binary mixing and its solution via the AMC is briefly 
reviewed. In Section 3, the Johnson-Edgeworth Translation is introduced with a highlight 
on the mathematical constraints associated with its application for the modeling of the 
mixing problem. Due to the previously established similarity of the JET frequencies 
with those based on the Pearson Family (PF), the Beta density of the first kind is 
also presented in this section. The PDF’s generated by these three models (AMC, 
JET and PF) are compared against each other and also with data generated by Direct 
Numerical Simulations (DNS) in Section 4. The results for the conditional expected 
scalar dissipation, and the conditional expected scalar dissipations for all the closures 
are discussed, respectively, in Section 5 and in Section 6. In Section 7, some theoretical 
remarks pertaining to the evolution of the scalar in an isotropic field are presented. With 
this presentation, the problems associated with all three closures become more clear. 
In Sections 2-7, the discussions are limited to those associated with the transport of a 
passive scalar from an initially symmetric binary state in isotopic turbulence. Therefore, 
in Section 8 some discussions are presented of the applications of the models for treating 
more general problems. This paper is drawn to a conclusion in Section 9. 

2 BINARY MIXING PROBLEM 

We consider the mixing of a scalar variable 4> = <£(x, t )(x is the position vector, and t 
denotes time) from an initially symmetric binary state within the bounds fa < <j> < <f> u . 
In this section, we assume that the lower and the upper bounds of the scalar range 
remain fixed ( Le . 4> u ,<f) e ll are constant). Within this domain, the single-point PDF of 
the variable 4> at initial time is given by 

Pi(4>, t = 0) = l -\ 6{<t> - 4> t ) + 6{4> - 4> u )l 


( 1 ) 
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FIGURE lb The comparison of the PDFs as predicted by the models with DNS data, (b) a 2 = 0.079. 


The closure problem in determining the PDF, Pi, is associated with the unknown 
conditional expected dissipation, e y and/or the conditional expected diffusion, Z>. These 
two are related through Eqs. (l)-(6), 


D(ht) 


1 gtfgi) 
Pi (<M) d<t> 


(7) 


At the single-point level neither the conditional mean dissipation nor the conditional mean 
diffusion are known (neither are their unconditional mean values). Their specifications 
require external information. 

With the application of the AMC, this external information is obtained in an implicit 
manner. As explained in detail by Pope (1991), the AMC involves a mapping of the 
random field of interest <j> to a stationary Gaussian reference field via a transformation 
<t> — O' Once this relation is established, the PDF of the random variable P\{<j>), 
is related to that of a Gaussian distribution. In a domain with fixed upper and lower 
bounds, Le. fixed <t>t, <f > u , the corresponding form of the mapping function is obtained by 
Pope (1991). The solution for a symmetric field with zero mean, < <f> >= 0 ,<f> u = 
is represented in terms of an unknown time r 


X(<foyT) = <t> u erf 



, G(r) = y/ exp(2r ) - 1. 


( 8 ) 


With this transformation, the PDF is determined from the physical requirement 


P\{x{<h,r),T)dx - PgOAoM^o, 


(9) 
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FIGURE k The comparison of the PDFs as predicted by the models with DNS data, (c) a 2 = 0.00247. 


where Pc denotes the PDF of a standardized Gaussian distribution, i.e. Pc{<M = 
-j-* txp(-<f>l/2 ). A combination of Eq. (8) and Eq. (9) yields 


Pi (far) 




( 10 ) 


In these equations, the relation between r and the physical time t is unknown in 
the context of a single-point description. This relation can be obtained only through 
knowledge of the higher order statistical properties of the scalar field. For example, it 
is shown by Madnia et al (1992); Frankel et al (1992a) that the mapping closure yields 
the algebraic relation for the normalized variance, 


< a 2 > (r) 

< a 2 > (0) 


2 

= - arctan 

7T 


( I 

VcVG 2 + 2 


)• 


( 11 ) 


in which the variance is related to the unknown mean dissipation e(l ) by integrating Eq. 

( 3 ), 

= _f(0 ’ (12) 


where, 








( 13 ) 
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3 JOHNSON-EDGEWORTH TRANSLATION 


The AMC captures some of the basic features of the binary mixing problem as described 
by Pope (1991). Namely, the inverse diffusion of the PDF in the composition domain 
from a double delta distribution to an asymptotic Gaussian distribution centered around 
< <f> >, as a 2 —* 0 (or G — * oc). This asymptotic Gaussian distribution near the mean 
scalar value cannot be realized in any of the previous mixing models based on the so 
called Coalescence/Dispersion (C/D) closures (Curl, 1963; Janicka et al , 1979; Pope, 
1982; Kosaly and Givi, 1987). And those modified C/D models which do yield such an 
asymptotic state, e.g. Pope (1982), do not predict the initial stages of mixing correctly 
(Kosaly, 1986). This deficiency of the C/D models in yielding asymptotic Gaussianity has 
been a motivating factor for recent investigations resulting in the development of the 
AMC (Pope, 1991). 

In the spirit of “mapping” to a specified reference field, it is speculated that there are 
perhaps other means of “driving” the PDF toward Gaussianity in a physically acceptable 
manner. In fact, this subject has been of major interest to statisticians and biometricians 
within the last century since the early work of Edgeworth (1907). The scheme was 
referred to by Edgeworth as the Method of Translation, and was later detailed by 
Johnson (1949a). In today’s literature of statistics and biometrics, the scheme is known 
as Johnson-Edgeworth frequency generation, and has many applications in statistical 
analysis. 

The essential element of Johnson-Edgeworth TYanslation (JET) is similar to that of 
the AMC. Namely, it involves the transformation of the random physical field, here, 0, 
to a fixed standard Gaussian reference field by means of a translation (or mapping) of 
the form 


<p = Z 


0o 

.7(0. 


7(f = 0) = 0 < 7(f) < 7(f — > oc) — ► oo. 


(14) 


In this equation, the function 7(f) plays a role similar to that of G in the AMC. With 
an appropriate form for the function Z, the scalar PDF is determined from Eq. (9). 
For application in the problem of mixing from an initially symmetric binary state of zero 
mean within a fixed domain 0/ = -0u < 0 < 0u, the appropriate JET must satisfy the 
following physical constraints: 


(i)Lim( 7 _»o)Z( — ) % H (0o). 
7 


(//)L/m (7 _ oo) Z(— ) % C0 O + O(0o) + 

7 

(iii) Z(^) is an odd function with respect to the scalar mean for any value of a 2 , (iv) 
Z(%) is bounded and is a non-decreasing function of 0 0, and — <f> u < Z < 0 tt at all 
times. 

(15) 

In these relations, H denotes the Heaviside function, and C is constant. Constraint (i) 
implies an initially symmetric and segregated binary state. The second constraint ensures 
an asymptotic Gaussian distribution for P\{<t>) near the mean scalar value. Condition (iii) 
preserves the symmetry of the PDF around the mean value at all times, and constraint 

(iv) implies the boundedness of the scalar field, Le. —<f> u < 0 < 0«- A function Z 
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which satisfies all the above constraints, is therefore expected to provide an acceptable 
means of approximating the PDF An example is the Logit - Normal (or tanh” ’-Normal) 
distribution, as originally proposed by Johnson (1949a). For the symmetric problem 
within a fixed domain, this distribution is produced by a mapping of the form 1 

Z = 4> u tanh (^) . (16) 

With this mapping, together with Eq. (9), the PDF of the scalar adopts the form, 


/M<M) 



2 r 


tanh 




(17) 


It is easily verified that this frequency satisfies the physical constraints of the symmetric 
binary mixing problem (Eq. (15)). At t =0, the PDF is approximately composed of two 
delta functions at 4> = ±<j> u > and as f cc the PDF adopts an approximate Gaussian- 
like distribution centered around the zero mean. These features are similar to those 
portrayed by the PDF generated by the AMC (Eq. 10)). 

This example demonstrates that with the satisfaction of the above indicated constraints, 
several other frequencies can be generated for effective modeling of the binary mixing 
problem. In fact, it is easy to show that the solution generated by the AMC can also be 
viewed as a member of the JET family. This is demonstrated by considering a translation 
of the form 



From Eq. (9), this translation yields the PDF 







(19) 


This frequency can be termed the erf” ’-Normal distribution and is identical to the form 
presented by Eq. (10). The difference is due to the terms containing G and 7 . But this is 
unimportant since in the context of single-point statistics neither of the two parameters 
can be determined by the PDF. Therefore, with G = both expressions are equivalent. 

With this equivalence, the closed form relation for the variance of the erf” 1 -Normal 
distribution has the same algebraic form given by Eq. (11). It is easy to show that many 
other distributions can be generated to display similar characteristics. In the discussions 
to follow, we only consider the Logit-Normal and the erf” 1 -Normal distributions, the 
latter being identical to the distribution generated by AMC. 


Pearson Family 

The similarity of the AMC and JET in generating equivalent PDF’s is also useful in 
explaining the applicability of the frequencies generated by the Pearson family (Pearson, 
1895). For a “bimodal” distribution, a physically acceptable frequency is the Pearson 


1 In recent literature, the Logit-Normal is usually expressed by the mapping Z - <p u 2[1 + expifoli)) 1 — 1. 
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FIGURE 2 Temporal evolution of the Logit-Normal PDF. 

Type I, known as the general form of the “Beta density of the first kind”. This density 
is typically expressed in a fixed domain within the range 0 < </> < 1, 

P.W = B J u ^-'(1 - < <f> < 1. (20) 

Here B denotes the Beta function, and the parameters j3\ and 02 are determined from 
the knowledge of the mean and the variance of the random variable. In a symmetric 
field within [0,1], < <\> >= j, 0\ = 02 = an ^ thus the PDF is characterized by the 
variance alone. 

The similarity of the Pearson distributions and the JET frequencies is well recognized 
in the statistics literature (see Johnson (1949a)). Therefore, with the equivalence of the 
AMC and the JET as demonstrated above, it is not surprising that the Beta density 
and the AMC are also similar. This similarity, without a mathematical proof, has been 
recognized in previous works (Madnia et al , 1991b; Madnia and Givi, 1992). 

4 COMPARATIVE ASSESSMENTS 

The probability distributions obtained from the three frequency generation methods 
described above are all capable of providing a reasonable stochastic approximation of 
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FIGURE 3a Temporal evolution of the centralized moments of the scalar variable as predicted by the 
models and the comparison with DNS data, (a) vs. t*. 

the mixing problem from an initially binary state. Namely, an approximate double delta 
distribution at / = 0, and an approximate Gaussian-like distribution as t — > oo. The 
former can be realized in the limit of unity normalized variance <r 2 (t)/a 2 ( 0) = 1. In a 
fixed composition domain, this corresponds to G = 0 for AMC (7 = 0 in JET), and 
to j3 - 0 for the Beta density. The latter is realized in the limit G, 7,/? — ► 00. The 
limiting Gaussian distribution for the AMC has been asserted by Pope (1991). For the 
JET, the criterion (ii) in Eq. (15) guarantees this condition. For the Beta density, the 
assertion of an asymptotic Gaussian distribution in the limit of zero variance is established 
in elementary texts on statistics (e.g. Casella and Berger (1990)). At the intermediate 
stages, however, the PDF’s are not identical. It is easily verified by Eqs. (10), (20) that 
the AMC and the Beta distributions become constant {P\{4>) — constant - \<j> u ) for 
G = (3 = 1. However, the Logit-Normal PDF does not yield a uniform distribution at 
any stage of its evolution. Also, as indicated by Johnson (1949a) it is not possible to 
provide a closed form algebraic expression similar to Eq. (11) for the variance of the 
Logit-Normal distribution. 

In order to make comparative assessments of the models, the frequencies generated 
by the three methods (AMC, JET, and PF) are compared with each other, and also with 
PDF’s generated by Direct Numerical Simulations (DNS). The DNS procedure is similar 
to that of previous simulations of this type. Since these simulations are not the major 
focus of this paper, only a brief outline of the procedure is described; for a detailed 
discussion we refer the reader to Madnia and Givi (1992). The subject of the DNS is 
a three-dimensional periodic homogeneous box flow carrying a passive scalar variable. 
The initial scalar field is composed of square waves with maximum and minimum values 
of 1 and 0, respectively. These limiting values are arbitrary, and can be translated to 
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FIGURE 3b Temporal evolution of the centralized moments of the scalar variable as predicted by the 
models and the comparison with DNS data, (b) /x$ vs. t*. 


appropriate 4> u , <f>t for comparison to each of the models. At time zero, in most regions 
of the flow, the scalar adopts these limiting values (with equal probability), and the 
spatial regions between the initial maxima and minima are smoothed by a polynomial fit. 
The degree of this smoothing is minimized, but limited by the computational resolution, 
to ensure an approximate initial double delta distribution. The computational scheme 
is based on a spectral-collocation procedure using Fourier basis functions developed 
by Erlebacher et al (1990a); Erlebacher et al (1987); Erlebacher et al (1990b). The 
hydrodynamic field is assumed isotropic and is initialized in a similar manner to that by 
Erlebacher et al (1990a); Passot and Pouquet (1987). The code is capable of simulating 
flows with different levels of compressibility (Hussaini et al , 1990). Here, only the 
results obtained for a low compressible case are discussed. The resolution consists of 
96 collocation points in each direction. Therefore, at each time step 96 3 is the sample 
size for statistical analysis. With this resolution, simulations with a Reynolds number 
(based on the Taylor microscale) of Re\ « 41 are attainable. The value of the molecular 
Schmidt number is set equal to unity. 

As indicated in Section 1, in order to compare the model predictions with DNS results 
a matching is required of the higher order statistics of the field as generated by each 
method. Here, this matching is done through the variance of the conserved scalar. These 

results are presented in Fig. 1. This figure indicates that at initial times, « 1, all 
the PDF’s are approximately composed of two delta functions at <j> - 0,1 indicating the 
initial binary state. At longer times, the PDF’s evolve through an inverse-like diffusion 
in the composition space. The heights of the delta functions decrease and the PDF’s 
are redistributed at other <j> values within the range [0,1], At very long times, the PDF’s 
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FIGURE 4a Comparisons of the normalized conditional dissipation as predicted by the three models with 
the DNS data, (a) <r 2 « 0.079. 


become asymptotically concentrated around the mean value in a manner that can be 
approximated by a Gaussian distribution. 

An interesting feature captured in Fig. 1(b), is the capability of the Logit-Normal 
distribution in depicting a subtle behavior in the frequency distribution. This feature is 
the double “hump” characteristic of the DNS data at intermediate times and cannot be 
realized by the AMC or the PF generated frequencies. All the previous DNS results 
including those of Eswaran and Pope (1988); Givi and McMurtry (1988); Pope (1991) 
portray this feature. The PDF’s generated by the AMC, and the Beta distributions 
adopt a constant value (of 1/2) when a 2 = ^ (for 0 < <f> < 1). This corresponds to 

G = 1,7 = n/2, (3 = 1. This uniform distribution is not exactly realized in any previous 
or present DNS results. Therefore, it can be speculated that in the absence of a better 
alternative, the Logit-Normal distribution may provide the simplest means of providing an 
assumed distribution for the statistical modeling of the symmetric binary mixing problem. 
The complete evolution of the Logit-Normal PDF is shown in Fig.2. 

Further quantification of the agreements noted above are made by comparing the higher 
moments of the scalar field. This comparison is made in Fig.3. In this figure, results are 
presented for the temporal variations of the kurtosis (/i4) and the superskewness (&,) 
of the scalar variable <f>. For the Beta density, the higher order moments are obtained 
analytically based on the knowledge of the variance. For the AMC, the analytical- 
numerical results by Jiang et al (1992) are used, while for the Logit-Normal PDF the 
moments are calculated strictly by numerical means. This figure shows that initially, all 
these moments are close to unity, and monotonically increase as mixing proceeds. For all 
the models, the magnitudes of the moments asymptotically approach the limiting values 
of 3 and 15, respectively, corresponding to those of a Gaussian distribution. The DNS 
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FIGURE 4b Comparisons of the normalized conditional dissipation as predicted by the three models with 
the DNS data, (b) a 2 = 0.013. 

results are in good agreement with the model predictions at all times. However, due to 
obvious numerical difficulties, the simulations could not be continued until the variance 
approaches zero identically. 

5 SCALAR DISSIPATION 

The results presented above indicate a good agreement between the model predicted 
single-point statistics (PDF’s and high-order moments) and the DNS data at all the 
stages of mixing. These results also suggest an approximate asymptotic Gaussian state 
for all the closure PDF’s and those of the DNS. Here, it will be demonstrated that the 
agreement between the DNS and the model predictions is very good at the initial and 
the intermediate stages of mixing. However, the agreement worsens at the final stages. 
Also it will be shown that none of the closures yield “exact” Gaussian distributions at the 
final stages of mixing. In doing so, it is useful to note that a Gaussian PDF is defined, 
and is only valid, for an unbounded domain. The frequencies generated here, are all 
defined within a fixed and finite domain. For AMC, it has been established (Gao, 1991; 
O’Brien and Jiang, 1991) that the finite boundary size at the initial time “maintains” 
its influence at all the subsequent stages of mixing. In other words, the PDF adopts a 
Gaussian distribution in the limit of zero variance only near the mean value of the scalar. 
In order to show the departure from Gaussianity at scalar values away from the mean, 
the conditional expected dissipation of the scalar field is considered. 

Given the PDF, as is the case here, Eq, (3) can be used to determine the expected 
conditional dissipation. It has been shown by Girimaji (1992) (and will be discussed in 
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FIGURE 5a Temporal evolution of the conditional expected diffusion normalized with the total dissipation, 
(a) erf - 1 ‘Normal. 

detail in Section 7) that for a valid PDF within a defined range, — <)> u < <f> < <t> u , the 
expected conditional scalar dissipation is given by 

'<*•'> = <?1 > 

where F denotes the cumulative distribution function (CDF) 

F(<M)= r Pi(rp,t)dip. (22) 

J —<t>u 

With Eq. (21), the expected conditional dissipation can be evaluated for a given PDF. For 
example, for a Gaussian distribution of zero mean, PG{<t>,° 2 ) = ex P(”§r> “°° = 
—a u < a < <r u = oo, with a non-stationary variance, a 2 = a 2 (/), it is easily shown that, 

\ + + w* £ s “ 

(23) 

Noting that <t> is an independent variable (of t), and evaluating the derivatives on the 
RHS of Eq. (23) yields, after some simple manipulations, 

. da 

£(<P,t) = constant = -a - — 


e(<M) = 


Pc(<t> 


4>,v 2 ) J dt \ 


( 24 ) 
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FIGURE 5b Temporal evolution of the conditional expected diffusion normalized with the total dissipation, 
(b) LMSE. 


with the implications (derived from Eqs. (12)-(13)), 


= constant = 1, (25) 

at all times for all phi values in the range -oc to +oo. Equations (24-25) indicate the 
independence of the conditional scalar dissipation and the composition domain for a 
Gaussian field. These results have also been obtained by Gao (1991); O’Brien and Jiang 
(1991) by following a different mathematical procedure. 

The conditional expected dissipation predicted by the models can be obtained by 
following a similar course. For the AMC and the PF distribution, the conditional 
dissipation fields have been obtained by Gao (1991); O’Brien and Jiang (1991) and by 
Girimaji (1992), respectively. For the purpose of the discussions to follow, these results 
are presented here in a different form for all three closures. For the erf" 1 -Normal 
distribution, the instantaneous CDF is given by 




1 + erf 


V2 


ert 


-'<£>])- 


(26) 


Therefore, with Eq. (21), the conditional dissipation can be expressed in terms of the 
corresponding PDF, 




Pt(+ 


1 d_ (If* 

4>j)dt \ 2 J _ 4 


erf 


~75 er f 'Or) 


Lv/2 


4 > u '\ 


dtp + -(< p + <Pu)\ 


(27) 
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Again, with an independence of <j> and t this equation reduces to, 


e(<j>,t) = - 


1 dH 
dt ’ 


15 


(28) 


where, 

For a PDF within a fixed domain, the integration procedure becomes simplified by 
evaluating the time derivatives inside the integral. In this way, the results can be expressed 
analytically. After some manipulations, 


dH 

dt 


jt(2 + 7 2 ) 


exp 


( 


-(1 + j ) 


erf 


-^f}' 


(30) 


and, therefore 


£(</>, 0 


tt7(2 + 7 2 ) 


exp 



erf 




(31) 


From this equation, the total dissipation is obtained by direct integration of the conditional 
mean dissipation field. The results, after significant algebraic manipulations yield 


e(0 = 


7r(2 4- 7 2 )\/4 + 7 2 


(32) 


e(<M) 


1 + sin 

2iW) 

‘(0 


1 — sin 

|53(0) J 


exp 


-2 


erf '( 


£>]'} 


(33) 


In the form presented above, Eqs. (31)-(33) portray several insightful features of the 
solution. First, Eq. (33) indicates that the conditional dissipation is always dependent 
on the magnitude of the scalar, and it maintains the same self-similar functional form 

of dependence exp 2 [ er / -1 (£o] J* This has been previously indicated by Gao 

(1991); O'Brien and Jiang (1991). Here, the amplitude e{<j> = 0, t) can be conveniently 
expressed in terms of the variance decay, which is very useful for further manipulations. 
Second, it is interesting to note the similarity of Eqs. (31) and (33) with the results 
obtained for the instantaneous dissipation of Fickian mixing of a conserved scalar in 
laminar non-homogeneous flows (such as the typical shear flows (Spalding, 1961; Linan, 
1974; Peters, 1984)). This similarity further asserts the “permanent” influence of the 
boundaries since in non-homogeneous mixing, the scalar bounds are “fixed” due to the 
physical constraints. Finally, Eq. (32) suggests an infinitely large dissipation at time zero, 
Le. when a 2 (t)fa 2 ( 0) = 1, and the asymptotic behavior 


Lim 


(<T : — >0) ~ I4>=0 - 1 


(34) 
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FIGURE 6a Comparisons of the normalized conditional diffusion as predicted by the three models and 
the LMSE closure, with the DNS data, (a) a 2 — 0.079. 


This limiting behavior near zero indicates the Gaussianity of the PDF only at the mean 
value of the scalar. 

Following a similar procedure, the conditional expected dissipation can be obtained 
for the other closures. For the Beta density in the range 0 < <f> < 1, the final results can 
be expressed as 

(35) 

where, / denotes the Incomplete Beta Function (Abramowitz and Stegun, 1972). For 
the Logit-Normal distribution, the corresponding form is 


0 


1 d 




(36) 


Neither of the equations (35-36) can be simplified further. Therefore, in order to evaluate 
the conditional expected dissipation (and the total dissipation), these equations must be 
evaluated numerically. 

In Fig. 4, the evolution of the conditional expected dissipation (normalized by the 
total dissipation) is presented for the models and the DNS data. This figure shows the 
similarity of the conditional expected dissipation for all of the models. The bell shape 
distribution is evident in all the figures with a maximum amplitude near the mean value. 
Also, as the variance decreases and the PDF becomes concentrated near the mean, the 
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FIGURE 6b Comparisons of the normalized conditional diffusion as predicted by the three models and 
the LMSE closure, with the DNS data, (b) a 1 = 0.013. 


amplitude tends to unity. This shape is typical of that observed in previous DNS results 
of Eswaran and Pope (1988); Nomura and Elgobashi (1992). 

The results in Fig.4 show the <j> dependency of the results at all the stages of mixing. That 
is, the PDF asymptotically adopts an apparent Gaussian-like distribution only near the 
mean value of the scalar, and the conditional dissipation does not become independent 
of the scalar everywhere. For the AMC, this has been discussed by Gao (1991); O’Brien 
and Jiang (1991). Considering the similarity of the three models, it is therefore concluded 
that all three models yield the same characteristics. These results also suggest a poor 
agreement between the model predicted conditional expected dissipations and the DNS 
data. Note that at the initial stages of mixing, the predicted results compare well with 
DNS data. However, with mixing progression, at smaller variance values, the agreement 
is only good near the mean scalar value and worsens near the bounds of the composition 
domain. This, as described above, is due to the permanent influence of the scalar 
boundaries at all the stages of mixing. 

6 SCALAR DIFFUSION 

Albeit directly related to the conditional expected dissipation (Eq. 7), it is useful to 
examine the behavior of the conditional expected diffusion in light of the discussions 
above. Given the PDF, again within the fixed range —<j> u < <j> < <f> u > the conditional 
expected diffusion can be determined from 


1 dF 
PiO?M) dt * 




( 37 ) 
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This equation is very useful in illustrating the properties of the PDF. For example, for a 
Gaussian distribution within an infinite domain 


dF 

dt 


and consequently 


<t> da 2 <t> 2 s 

VW dt CXP( ty )% 


£>(<M) _ 

e{t) <T 2 (t)' 


(38) 


(39) 


It is noted that Eq. (39) is in accord with the Linear Mean Square Estimation (LMSE) 
closure (O’Brien, 1980). 

The mean conditional diffusion can be determined for the three models considered. 
For the erf -1 -Normal PDF with zero mean 


dF 

dt 


1 d^i 
y/2n dt 




* r f & 

<Pu 


(40) 


Again with explicit equations for the total dissipation and the variance, it is possible to 
obtain an algebraic expression for the conditional expected diffusion. The results after 
substantial algebraic manipulations yield 



t) 

I -sfr 

1 + sin 

™ 2 (0 
L 2<x'(0) 


<0 

\ sin [i£$] \ 

1 - sin 

[ 2<r*(0) 


exp 






(41) 

In this form Eq. (41) is very pleasing since it does reveal the (t,tf>) separability, and thus 
the self-similarity, of the diffusion field. The terms inside the parenthesis on the RHS are 
time dependent, whereas the remaining terms depend explicitly on <f> only. As indicated 
by O’Brien and Jiang (1991), this separability cannot be easily deduced from Eq. (5), 
but is possible with the analytical procedure followed above. The temporal evolution of 
the conditional expected diffusion for the erf _1 -Normal distribution, and its comparison 
with that of the LMSE closure is presented in Fig.5. 

By following the procedure above, analogous expressions are obtained for the other 
two closures. Namely, 


g(jM) _ 2<f>u d^t 

<t) 



<t> 

tanh *(-“), 
<Pu 


(42) 


for the Logit-Normal distribution of zero mean, and 

g(^0 = 2 a/,CM) 

€(0 PMj) da 2 


(43) 


for the symmetric Beta density within [0,1]. Equations (42)-(43) cannot be simplified 
further due to the lack of an explicit analytical relation for the variance of the Logit- 
Normal distribution (Johnson, 1949a), and the unknown analytical form of the derivatives 
of the Incomplete Beta Function. 
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FIGURE 7 The temporal variations of <£m a x and ^ mjn generated by DNS. 

In Fig. 6, results are presented of the conditional expected diffusion as predicted by the 
three models, and also that of the LMSE closure. In these figures, the DNS data are also 
provided at several variance values. The similarity of the modelled results are once again 
revealed in these figures, which is expected in view of the PDF similarities. At all times, 
the conditional diffusion field has an odd distribution near the mean scalar value. On 
the right half of the composition domain, all three closures yield a monotonic decrease 
of D to an instantaneous minimum, and then a monotonic increase to zero at the upper 
bound of the scalar. The location and the magnitude of the instantaneous maxima and 
minima is not the same for the three closures. Also, as Eqs. (41)-(43) indicate, the 
zeroes of D can only be realized at <f> = 0, ±<f > u . At the initial times, Le. large variances, 
all three closures agree reasonably well with the DNS data. This agreement is better 
for the three models than for the LMSE closure. However, as the variance becomes 
smaller, the agreement between the model predictions and the DNS data worsens. It 
is noted that as the variance becomes small, all the closures yield a Gaussian-like PDF 
near the mean value of the scalar. This is shown in the figures near (j> =< <t> > (= j 
for DNS), where the predicted results are in accord with the LMSE closure, Le. linear 
profiles of similar slopes. In this region, the results are also in accord with DNS data for 
all the closures including the LMSE. However, again, the predicted results deviate from 
the DNS data away from the mean value. It is clearly noted that the DNS generated D 
values do not go to zero at the scalar bounds. 

7 EVOLUTION OF THE SCALAR HELD 

The problems described at the conclusion of the previous two sections stem from a 
lack of capability of all of the models in accounting for the variations of the scalar 
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FIGURE 8 The Beta density (Pearson Type II) for a domain with moving boundaries, and /? — 0.1. 

bounds as the mixing proceeds. For all three models, the PDF is always defined within a 
fixed range through its course of evolution. It is easy to show that both the conditional 
expected dissipation and the conditional expected diffusion are correctly predicted by all 
the models near the mean scalar value. For the erf - ^Normal distribution, this is evident 
from Eqs. (33) and (41) and can be also shown by analyzing the behavior of Eq. (31) 
near the region <f> % 0, as the variance becomes small. Noting that 



Li m^ 0) er fA a 

<p u L <Pu 

(44) 

and, from Eq. (11) 

dl —2<t>u da 

(45) 


Lim W^o) dt ~ 

it is easily concluded that 




Lim (t ^ oo)f(4> % 0,0 = 

dt 


(46) 
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Following the same procedure, it is derived 

Lim {t ^ oo) D{4> ss 0,/) = — (47) 

(j* 

Due to the similarity of the three closures, it is reasonable to expect similar behaviors 
for the other models as well. Equations (46)-(47) indicate a Gaussian-like distribution 
near the mean <f> as 0 (Eqs. (24,39)). This is in accord with the DNS data. However, at 
distances away from the mean value the predicted results do not correspond to that of 
a Gaussian field. Neither do these results agree well with DNS data. The deficiency of 
the models in predicting the DNS results is made clear by considering the bounds of the 
scalar field as the mixing proceeds. This is demonstrated in Fig. 7, showing the temporal 
decay/growth of the scalar maxima/minima obtained by DNS. This trend is consistent 
with physical intuition, but is not incorporated into any of the three models. In the 
AMC and the JET generated frequencies, due to the nature of the translation Z{4> o, t) 
and the constraints imposed by Eq. (15), the scalar is always bounded within the same 
range. This problem is also encountered in the PF, in that Type I and II distribution 
families are always defined within the same fixed domain regardless of the magnitude of 
the variance. 

With the examination of the PDF transport equation, it is shown that the physics 
of the problem requires the migration of the scalar bounds toward the mean value 
as the mixing proceeds. That is, the instantaneous values of the scalar minima and 
maxima change with mixing progression. To demonstrate this, again consider a symmetric 
field with a PDF, P\(<j>,t), defined within the time-dependent domain of zero mean, 
<p£[<t>min(t ) = -<£max(f )< 0max(O] At all the stages of the evolution, the PDF must satisfy 
the physical requirements 

r^max (f ) 

/ P\((j>,t)d<j> ss 1, 

d <Amax(f ) 

r<t>mu0) 

<<P>= / <j>Pi(4>,t)d<t> = 0, 


a\t) = 

J<p„ 




The first of Eq. (48) requires 


d_ r* 
dt L 


P\ (4>, / )d (p = 0 


Evaluating this integral via Leibnitz’s rule, and making use of Eqs. (3)-(7), it is shown 
that 


p.(*_( 0,0*37= ’ - {£[.(*,<> i>.(* 
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Following the same procedure for the second of Eq. (48), yields the obvious requirement 


r4> m*x(0 

/ D(<j>,t)P\(<f>,t)d4> = o. 




The third part of Eq. (48) yields Eq. (12), and 


(51) 


PM max (0 max(t ), t ) — P i(<^min(t ), t )f(^rnin(t ), t ) — 0, 


PM max (0.0 
^ 1 (<^min(f )»0 


d 4>mzx 

dt 


^(^max(f )> 



dt 


)i 0 


= 0 . 


(52) 


The remaining parts of Eq. (48) yield higher order statistical information pertaining 
to the inner integrated evolution of the conditional expected dissipation and diffusion, 
and their relation with the higher central moments. With an additional assumption of a 
nonzero PDF within the region of its definition, that is by defining 0 m ax(O and <£min(0 
as the extreme locations with nonzero PDF, a combination of Eqs. (50) and (52) yields 


£(0max(f )» 0 — £(0min(f )i0“®i 

- D(A,„(0,'), 

0max(O) = 4>u i 0min(O) = <f>(. (53) 

Equation (53) indicates that with fixed boundaries, the conditional dissipation would 
adopt a zero slope at the boundaries and the conditional diffusion would also be zero 
there. However, Fig. 7 indicates that in a physical situation the boundaries are not fixed 
and move inwards as the mixing proceeds. It is interesting to note that this problem 
is not observed in the numerical results obtained by the C/D type closures. That is, 
while the C/D closures are not capable of predicting the PDF evolution in accord with 
DNS data, they do have the mechanism for shrinking the bounds of the composition 
space. Obviously, in the context of single -point description without the knowledge of 
the dissipation field, it is not possible to determine a priori the temporal bounds of the 
scalar field. Therefore, the closures can be modified only by making further assumptions 
in describing this transport. For a general case, the JET frequencies can be generated 
by the original form proposed by Johnson (1949a) 

<X<t>o,t) = H‘)z +0(0, (54) 

where the additional parameters A(r) and g(t) provide the extra degrees of freedom in 
order to account for the variations of the instantaneous boundaries of the composition 
domain. For the PF, the problem can be overcome, for example, by considering a 
“four-parameter Beta distribution” 
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FIGURE 9a The comparison of the conditional expected scalar dissipation normalized with the total 
dissipation with DNS data as predicted by the AMC with the scalar bounds determined from the DNS results, 
(a) a 2 = 0.079. 
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0 0min(f ) 
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i $min {§ ) — 4* ^ 4* max (0, 


(55) 


with the extra two parameters being <£ m m(0 and </> m ax(0* F° r a symmetric PDF in the 
range [0,1] = [</> m j n (/ = 0) = </>/, </> ma x(f = 0) = <f> u \ therefore, the variance decay 
can be influenced by increasing 0, and/or by decreasing the scalar range A 4>{t) = 
0max(O — 0min(O* The former recovers the well-known two-parameter Beta distribution 
(Pearson Type II), while the latter is approximately equivalent to the LMSE closure 
(O’Brien, 1980). This latter case is presented in Fig. 8 showing a symmetric Beta density 
with 0(t) = fixed = 0.1. Note that as the mixing proceeds, the variance decays but the 
PDF preserves its initial approximate double delta shape. In a physical problem, the 
situation is somewhere between these two limiting cases. The exact situation depends 
on the characteristics of a particular mixing problem. 

The discussions above suggest that in order to predict the final stage of mixing correctly, 
the effects of mixing on the shrinkage of the domain must be taken into account. To 
demonstrate this point, the results shown in Fig. 7 can be incorporated into the mixing 
models to determine the evolution of the conditional expected dissipation and diffusion. 
This is done here only for the erf" ’-Normal, and the results of the conditional expected 
dissipation are shown in Fig. 9. In the calculations resulting in this figure, analytical 
solutions are not possible for the moving boundary case. This is demonstrated by the 
equivalent form of Eq. (29) 
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FIGURE 9b The comparison of the conditional expected scalar dissipation normalized with the total 
dissipation with DNS data as predicted by the AMC with the scalar bounds determined from the DNS results, 
(b) <r 2 = 0.013. 
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With this equation, therefore, the effect of the temporal variation of the PDF on the 
conditional dissipation is through the dH jdt term in Eq. (28). This term has the form 
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The first term on RHS of Eq. (57) is the same as that in Eq. (30), and the effects of moving 
boundaries manifest themselves through the second term. This term cannot be evaluated 
analytically. However, Eqs. (56)-(57) show that due to the direct dependence on 
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the conditional dissipation does not retain its original functional dependence, suggested 
by Eq. (33). Also, Fig. 9 shows that the effect of the moving boundaries is to force the 
conditional expected dissipation to zero at the current scalar maxima/minima. Therefore, 
the predicted results compare much better with DNS data than those presented in Fig. 
4. Due to the similarity of the PDF’s, it is expected that the other two closures would 
also behave in the similar fashion. 

The influence of boundary encroachment is also sensed in the conditional expected 
diffusion field. For the erf -Normal scalar PDF, the equivalent form of Eq. (41) is 


^ ,, 2 <t> max (t)d'r 


erf ’ Gm«(o) + (ui) 
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with an average dissipation of 
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Equations (58)-(59) show the influence of the boundary movement through the last term 
on the RHS of both these equations. With these additional terms, the normalized form 
similar to Eq. (41) is not very useful, and Eqs. (58)-(59) are evaluated numerically. 
The equivalent of Eq. (58) for the Logit-Normal and the Beta density are, respectively, 




tanh- 1 (— £— ) + (— £— ) 
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(61 > 

An interesting characteristic displayed by Eqs. (58) and (60) is the influence of the 
terms containing the temporal derivative of ^ ma x(0* N° te that at the boundaries, Le. 
$ = 4> max(0» term on l ^ e RHS of these equations vanish, but the last term 

prohibit the conditional expected diffusion from going to zero. This is in accordance 
with the DNS data as shown in Fig. 6. In order to demonstrate this more clearly, results 
are presented in Fig. 10 of the conditional expected diffusion predicted by the erf -1 - 
Normal, with the input of the variance and the scalar bounds from DNS. A comparison 
between this figure and Fig. 6 show the influence of the boundary movement, and a 
better agreement between the model predictions and the DNS data. This agreement 
is more pronounced at scalar values away from the mean. Near the mean value, the 
influence of the boundary migration is slight, as also indicated by Eqs. (58) and (59). 


8 DISCUSSIONS AND APPLICATIONS 

In the previous section, a rather detailed discussion was presented of the problem of 
scalar mixing from an initially symmetric binary state. These discussions were primarily 
intended to provide a means of assessing the differences between the currently available 
tools in probability modeling of the scalar mixing problem. This problem is of significant 
interest, considering the extent of previous works focused on its analysis (Pope, 1976; 



i->' j ruLliifm>9. iNt. 

TeX output 1993 01 21:1426 (26) 


26 


R. S. MILLER, S. H FRANKEL, C. K. MADNIA AND P. GIVI 


10.0 



FIGURE 10a The comparison of the conditional expected scalar diffusion normalized with the total 
dissipation with DNS data as predicted by the AMC and the LMSE closures with the scalar bounds 
determined from the DNS results, (a) a 2 - 0.079. 

Pope, 1979; Pope, 1982; O’Brien, 1980; Dopazo, 1973; Kosaly and Givi, 1987; Pope, 1991; 
Gao, 1991; O’Brien and Jiang, 1991; Nomura and Elgobashi, 1992). The results obtained 
here are particularly useful in highlighting some of the deficiencies of these closures, 
and in suggesting future research towards overcoming these drawbacks. There are, 
however, many other physical problems that are not subject to the restricting conditions 
imposed in these analyses. In this section, therefore, some discussions are presented as 
to the practical implications of these models, together with some speculations on their 
extensions for future applications. 

Perhaps one of the most important practical applications of the closures considered here 
is the treatment of reactive flow phenomena. In fact, the most important advantage of 
scalar PDF methods is due to their applicability in the modeling of turbulent combustion 
(Pope, 1979; Pope, 1985; Pope, 1990; Kollmann, 1990; O’Brien, 1980). The results 
generated here can be used directly in the modeling of mixing controlled homogeneous 
chemically reacting systems. Namely, in examining the compositional structure of a 
reacting system under chemical equilibrium, or in determining the limiting rate of 
reactant conversion in a simple chemistry of the prototype Fuel -f Air — ► Products. 
The determination of this rate has been the subject of extensive investigations over 
the past forty years (see Hawthorne et al (1949); Toor (1962); Williams (1985)). It is 
now well-established that in a mixing controlled binary irreversible reaction of this type, 
the statistics of the reacting fields can be related to those of an appropriately defined 
conserved scalar (such as <j > ) (Bilge r, 1980; Toor, 1975; Williams, 1985). Therefore, 
the frequencies generated herein can be utilized for estimating the statistics of the 
reacting field with an infinitely fast chemistry model in a homogeneous flow with an 
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initially segregated reactants under stoichiometric conditions. Albeit very restricting, this 
problem is of great practical importance for modeling and design of batch mixers and 
plug flow reactors in which these conditions prevail (Toor, 1975; Brodkey, 1981; Dutta 
and Tarbell, 1989). Madnia et al (1991b); Madnia et al (1992) have shown that with 
the erf^-Normal (AMC) and the Beta density models, this rate can be predicted by 
simple analytical means. For the Logit-Normal density, a complete analytical solution 
cannot be obtained and determination of the statistics requires numerical integration of 
the PDF. The results generated by these closures agree with DNS data better than those 
obtained by means of the C/D closures (McMurtry and Givi, 1989), or other models 
previously available in the chemical engineering literature (Dutta and Tarbell, 1989) (see 
Givi (1989) for a review). Also, the results provided by the AMC (Frankel et al, 1992a) 
are shown to compare well with experimental data on plug flow reactors if the additional 
information pertaining to the evolution of the scalar length scale is accurately provided. 

The most obvious issues in regard to the applications of these models are associated 
with their extension for the treatment of (1) non-symmetric binary scalar mixing, (2) non- 
binary scalar mixing, (3) multiple scalar mixing, and (4) non-homogeneous mixing. The 
first problem constitutes a more general form of the binary mixing problem and is also 
important for the analysis of non-stoichiometric reacting systems. The second problem is 
appropriate for the analysis of other mixing systems in which the initial conditions are not 
of a two-feed configuration. The third problem is of interest in reacting systems in which 
the transport of a passive scalar (like <f>) is not sufficient for a complete analysis. For 
example, any reacting system under non-equilibrium conditions. Finally, the importance 
of the fourth problem is obvious in view of the fact that the flow within most practical 
mixing devices cannot be assumed homogeneous. 

In regard to the first issue, all of the three closures considered here can be used for the 
probability modeling of scalar mixing within a fixed scalar domain. The use of the AMC 
is straightforward, but the mathematical procedure is somewhat complex (Madnia et al, 
1992). The Pearson frequencies can be generated easily for non-symmetric problems. 
In this case, the Pearson Type I provides a reasonably accurate representation of the 
scalar field regardless of the degree of asymmetry of the PDF (Frankel et al, 1992b; 
Madnia et al, 1992). The use of the JET in this regard is most difficult, since closed 
form analytical expressions are not available for the variance of the scalar by which the 
PDF can be characterized (Johnson, 1949a). In treating these problems, therefore, the 
first two models can be more readily employed and subsequently used for the treatment 
of mixing controlled reacting systems under non-stoichiometric conditions. In fact, as 
demonstrated by Madnia et al (1992) the solution of the non-symmetric form of the 
AMC and the Beta density provide a very good means of predicting the limiting rate 
of reactant conversion in homogeneous reacting flows. However, it should be indicated 
that with both models the problem associated with the scalar bounds still exists and must 
be dealt with as discussed in Section 7. 

In addressing the second issue, it is obvious that the AMC is more appropriate than 
the other closures for simulating the mixing problem from an initially “arbitrary” state. 
The extension of JET and PF for treating multi- (higher than bi-) modal distributions 
have been reported in statistics literature. However, as the degree of modality of the 
PDF increases the procedure becomes more complex and not suitable for practical 
applications. Fortunately, in most mixing problems in simple flows, Le. homogeneous 
turbulence and turbulent shear flows, the PDF exhibits strong bimodal features (Madnia 
et al, 1992; Frankel et al, 1992b). In those cases, the use of the Beta density can be 
justified. In fact, in non-homogeneous flows it is easier to use this density, at least until 
further developments of the AMC for practical applications (see Frankel et al (1992b); 
Gaffney et al (1992)). 
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FIGURE 10b The comparison of the conditional expected scalar diffusion normalized with the total 
dissipation with DNS data as predicted by the AMC and the LMSE closures with the scalar bounds 
determined from the DNS results, (b) a 2 = 0.013. 


The extension of all of the three models in describing multi-scalar mixing is possible. 
The problem naturally falls within the realm of the multivariate statistical analyses. In these 
analyses, the implementation of the AMC is relatively straightforward since it provides 
a transport equation for the joint PDF’s of the scalar variable in a multivariable sense 
(Pope, 1991). However, it is not presently clear how to devise an efficient computational 
procedure, typically based on Monte Carlo methods (Pope, 1981), for the numerical 
treatment of these equations. Some work in this regard is currently under way (Valino 
and Gao, 1991). The extension of assumed distributions based on the Beta density for 
treating multi-scalars is more straightforward but less trivial to justify. The most obvious 
means is to implement the multivariate form of the PF. The direct analog of the Beta 
density is the Dirichlet frequency (Johnson, 1987; Narumi, 1923; Johnson and Kotz, 
1972). The application of this density in modeling of multiple species reactions has been 
discussed by Girimaji (1991a); Girimaji (1991b); Gaffney et al (1992). However, the 
use of the Dirichlet frequency cannot be justified for modeling of reacting flows in a 
general sense (Frankel, 1992). Finally the extension of the JET in generating multivariate 
frequencies has been reported in statistics literature since the subsequent work of Johnson 
(1949b). As one may suspect, the procedure is more complex, and the same reservations 
as those associated with the Dirichlet distributions apply. 

All of the models considered here can be extended for the analysis of non-homogeneous 
mixing (and reacting) systems. Obviously, in most cases, the problem requires numerical 
integration of the appropriate conservation equations. For instance, the AMC can be 
invoked in the mixing step of a fractional stepping procedure, similar to that of typical 
Monte Carlo procedures (Pope, 1981). The PF densities (e.g. Beta or Dirichlet) and JET 
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generated frequencies require modelled transport equations for the first two moments 
and cross moments of the scalar field. These equations, “hopefully”, include the essential 
information pertaining to the spatial inhomogeneity of the flow. Naturally, the PDF is 
not generally symmetric, and must be determined from the knowledge of the parameters 
/?!,/? 2 , 7? A, q , and the local </> max (f ), values. With this information, all the higher 
order statistics of the scalar field can be determined. In regard to this last issue, it must 
be indicated that the Beta density has been extensively used for the modeling of non- 
homogeneous reacting systems (e.g. Rhodes (1975); Jones and Priddin (1978); Lockwood 
and Moneib (1980); Janicka and Peters (1982); Peters (1984); Frankel et al (1992b); 
Gaffney et al (1992); for recent reviews, see Givi (1989); Priddin (1991)). Due to their 
special mathematical properties, the Beta and/or the Dirichlet frequencies yield relatively 
simple analytical solutions for the higher order statistics of the reacting fields. From 
this point of view, the use of the PF is more practical than the AMC since the solution 
procedure does not require the numerical treatment of the PDF transport equation. This 
point has been discussed in detail by Girimaji (1991b). However, as indicated above, the 
use of the Dirichlet frequency cannot be justified for modeling of unpremixed reacting 
flow in a general sense. Also, there is no way of implementing this density directly for 
modeling of non-equilibrium flames, involving strong correlation of the temperature and 
the species mass fractions. Even with the assumption of statistical independence of the 
reacting species and the temperature, the question of the local scalar range imposes a 
severe restriction on the validity of this approximation. For example, it is demonstrated 
by Gaffney et al (1992) that in the modeling of a reacting turbulent shear flow, depending 
on the a priori choice of the magnitudes of the local scalar bounds the predicted results 
can be altered significantly. Obviously, this problem is not eliminated with the usage of 
JET frequencies in either a univariate or multivariate sense. 

9 CONCLUDING REMARKS 

It is shown that the family of frequencies generated by the Johnson-Edgeworth Translation 
(JET) provides a reasonable means for statistical modeling of binary symmetric scalar 
mixing in homogeneous turbulence. It is also shown that the results predicted by one 
of the members of this family is identical to the solution generated by the Amplitude 
Mapping Closure (AMC) of Kraichnan. This similarity is useful in two regards: (1) 
establishing a mathematical reasoning for the similarity of the probability frequency of 
the Pearson Family (PF) and that of the AMC for the description of the problem, and (2) 
suggesting the possible use of other members of the JET frequencies in approaches in 
which the Probability Density Function (PDF) is assumed a priori. The PDF’s generated 
by all these models are shown to compare well with each other and also with the results 
obtained by Direct Numerical Simulations (DNS). However, none of the models are 
capable of accurately predicting the conditional expected dissipation and the conditional 
expected diffusion of the scalar field. This problem is associated with the incapability of 
the models to account for the migration of the scalar bounds as mixing proceeds. 
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4 Appendix II 


This work summarizes all our efforts to-date in LES of reacting turbulent flows. The primary 
conclusion drawn from this work is that the procedure used in k transport equation should 
be also followed for the g-transport equation. We are presently in the process of undertaking 
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Abstract 


A priori and a posteriori analyses of homogeneous and non-homogeneous flows are conducted 
in order to assess the validity of assumed Probability Density Function (PDF) models as po- 
tential subgrid scale (SGS) closures for Large Eddy Simulations (LES) of turbulent reacting 
flows. Specifically, a priori analysis of homogeneous and non-homogeneous turbulent flow 
is conducted, for both equilibrium and non-equilibrium chemistry, to investigate the po- 
tential of the Pearson Family PDF’s (Beta and Dirichlet) as SGS models. Also LES of a 
two-dimensional turbulent reacting shear layer are conducted using a hybrid one-equation 
Smagorinsky/PDF SGS closure model. The traditional Smagorinsky closure, augmented by 
the solution of the subgrid turbulent kinetic energy (TKE) equation, is employed to account 
for the hydrodynamic fluctuations and an assumed PDF is employed for treatment of the 
scalar fluctuations. An isothermal reaction of the type A + B — > Products is considered. For 
this reaction, the assumed PDF approach requires the local mean scalar values and the local 
subgrid covariance. This covariance is obtained by solving an additional modeled transport 
equation. This approach results in simple algebraic closures for the chemical source terms 
appearing in both the species continuity equations as well as in the species covariance equa- 
tion. The simulations are performed by using a hybrid spectral/finite-difference numerical 
algorithm. Results are compared to those obtained via Direct Numerical Simulations (DNS) 
of the same flow in order to assess the validity of our hybrid S/PDF closure. 
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1 Introduction 


Large Eddy Simulations (LES) lie somewhere between Direct Numerical Simulations (DNS) 
and Reynolds Averaged Navier-Stokes (RANS) computations. In DNS all relevant scales 
are resolved and no attempt is made at modeling the statistical behavior of the flowfield 
(Givi, 1989). In RANS calculations, all scales are modeled and moments above a certain 
level (usually second) are truncated (Launder and Spalding, 1972). In LES the large, energy- 
containing scales of motion are treated directly, while the effect of the small scales is modeled. 
Thus scales below a certain level are truncated and modeled. LES is currently viewed as a 
research tool, but because of the savings in computational resources over DNS, it has the 
potential to be used for engineering applications. 

Over the past 30 years, since the early work of Smagorinsky (1963), there has been rela- 
tively little effort, compared to that received for Reynolds Averaged Navier-Stokes (RANS) 
calculations, to make full use of the LES approach for engineering applications. The most 
prominent model has been the Smagorinsky eddy viscosity model which relates the unknown 
subgrid scale Reynolds stresses to the large scale rate of strain. The eddy viscosity performs 
the predominant role of mimicking the dissipative behavior of the unresolved small scales. 
One-equation models, which provide the subgrid scale turbulent kinetic energy as the veloc- 
ity scale for the eddy viscosity, have shown improvement over algebraic models for coarse 
grid simulations and for modeling higher order moments (Lilly, 1967; Schumann, 1975; Ho- 
riuti, 1985). In transitional flows where the assumption of a balance between production 
and dissipation may not always be valid, higher order models allow more freedom for the 
subgrid scale eddies to adjust to local flow conditions (Rogallo and Moin, 1984). 

LES of mixing layers have been conducted by Mansour et a I. (1978), Cain et a/. (1981), 
Lesieur et a/. (1988), Maruyama (1988), and more recently by Ragab et a/. (1992). With 
regards to applying LES to chemically reacting turbulent flows the literature is even more 
scarce. Schumann (1989) seems to have been one of the first to conduct LES of reacting 
flows. However, his assumption to simply neglect the subgrid scale (SGS) contributions from 
the chemical reactions is questionable. The importance of the scalar fluctuation correlations 
in turbulent reactive flows has been well recognized (Libby and Williams, 1980). The success 
of PDF methods for the treatment of chemically reacting flows has been associated with the 
closed form of the chemical source term (Pope, 1985; 1990; Givi, 1989). These methods 
have been used extensively for the modeling of turbulent reacting shear flows in a variety ol 
different RANS applications (Lockwood and Naguib, 1975; Bockhorn, 1987; Frankel et a/., 
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1990). Recently it has been demonstrated that the Beta density provides reasonable means 
for predicting the statistical behavior of a conserved scalar in both homogeneous and non- 
homogeneous flows (Madnia and Givi, 1992; Madnia et a/., 1991; Frankel et al, 1991,1992). 
Girimaji (1991) suggests the use of a Beta density for modeling of turbulence/chemistry 
interactions in multivariate species fields. Gaffney et a I. (1992) and Baurle et al. (1992) 
have used the joint Beta PDF in RANS calculations with complex kinetic mechanisms. The 
idea to hybridize the DNS approach with PDF methods for treatment of reacting flows is a 
natural extension of the current technique and is warranted in view of predicted computer 
capabilities for the foreseeable future. 

In this work, we first consider a priori analysis of a homogeneous reacting flow in order to 
assess the extent of validity of the Pearson type PDF’s as a SGS model. Then, we propose 
a hybrid SGS model in order to conduct LES of turbulent reacting flows. In both cases, 
a chemical reaction of the type A + B — > Products is considered. The model consists 
of the traditional Smagorinsky closure for the hydrodynamics with the turbulent kinetic 
energy providing the velocity scale and an assumed PDF approach for treatment of the 
chemical reaction rate terms. This technique requires the solution of an additional transport 
equation for the subgrid scale turbulent kinetic energy and the species covariance. Results 
are compared to those obtained without the PDF chemistry model as well as those obtained 
directly by DNS. 


2 Hydrodynamic Formulation 


The incompressible Navier-Stokes equations are, 



(i) 


du{ d(uiUj) 1 dp d 2 U{ , oX 

1 — = h v - — ( l ) 

dt dxj pdx t dxjdxj 

Typically one defines a filter in order to delineate the resolved or large scale field from the 
subgrid scale motions, 

u = J G ( r — r')u(r') dr' ( 3 ) 

where G is a filter function with characteristic length A. There are two main approaches 
towards this filtering process. They are the prefiltering approach and the grid averaging 
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approach . In this work we choose to employ the latter for its ease of implementation. Thus 
all flow variables are considered to be grid averaged variables and as such can be decomposed 
into subgrid scale fluctuations from the large scale field, 

u'i = Ui-TTi ( 4 ) 


Applying the filter to the governing equations, 

du 7 


dxi 


dui d(uiUj ) 
dt + d Xj 


= 0 


1 dp d 2 U{ 
+ v- 


pdxi dx } dxj 


The nonlinear term, looks like 


V~Uj = (u'i + U7)(u' + Uj) = u'iUj + u'Uj + U.U' + Ui Uj 


( 5 ) 

( 6 ) 

( 7 ) 


The last term depends on large scale components and is computable in LES. The subgrid 
scale Reynolds’ stresses are defined as, 


R,j ~ u\u' } + u[uj + iiiu'j 


( 8 ) 


This forms the focus of hydrodynamic SGS modeling. Typically one decomposes the SGS 
stress into the sum of a trace-free tensor and a diagonal tensor, 


Rij = (Rij — -SijRkk) + ~$x]Rkk = Tij -f -SijRkk 


( 9 ) 


Substitute this into the equations, 


or defining, 


du{ d(ui uj) _ 1 dp drij d(\8 l3 Rkk ) ^ <9 2 u, 

dt dxj pdx{ dxj dxj dxjdxj 


— pi 

P = -+ -Rkk 
P 3 


( 10 ) 


( 11 ) 


and employing the Smagorinsky eddy viscosity model for closure of the Reynolds’ stress, we 
have: _ . 

n , ; = 2^S is = n (12) 


4 



In the original Smagorinsky model the eddy viscosity, i/t, was related to the large scale 
strain rate. A suggestion made by Kwak et al. (1979) and discussed by Mansour (1981) is to 
relate the eddy viscosity to the trace of the resolved vorticity field. This is done in the hope 
of dealing with turbulent/non-turbulent interfaces more effectively, because the vorticity is 
zero in regions of irrotational flow, whereas the strain rate may not be. This model, while 
adequate for closure of the first order large scale transport equations, was found insufficient to 
achieve closure of the second order SGS species covariance equation. Therefore, the decision 
was made to go to a one-equation hydrodynamics closure model, specifically we solved the 
transport equation for the SGS turbulent kinetic energy, k = |(u'u;). The eddy viscosity is 
of the form, 

v r = C k ^\fk (13) 


where C k is a model constant chosen as 0.01, and A is the filter width chosen so as to enable 
comparison to the same physical space as in the DNS calculations. Then our modeled 
momentum equation becomes, 


du % 

dt 


+ 


dim u j ) 

dxj 


+ 4 . dv T duj 

pdxi dxj y T dxjJ dx : dx t dxj 


(H) 


The form of the modeled TKE equation is shown below, 


Dk 

~Dt 


d_ 

dx t 


(v 4 ) o 

V (Tk OXi , 


jfcf 

+ P k ~ C D — 


(15) 


where P k is the production term given by, 


P k = Vj 


(duj dui 
dxi dxj 


du] 

dxi 


(16) 


This term can be rewritten in terms of the large scale strain rate, Sij = | ijfal 4" If)’)’ 
as P k — 2urSijSij. The form of this term arises from extending Bousinesque’s relation 
for the laminar stress tensor to the turbulent stress tensor. This relation postulates that 
the the unknown subgrid stress tensor is proportional to the large scale strain rate, the 
proportionality being an eddy viscosity. We found that this traditional modeling resulted in 
too small a production of TKE at the mid- and far- field regions of the flow domain. This 
will be shown in the results section, but this problem necessitated an unorthodox modeling 
approach for the production term in the TKE. Thus, we found out that if we modeled the 
production term using the square of the rotation tensor we get a much improved solution 
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for the TKE over more than two thirds of the flow domain. This model is based on two 
observations: (1) the regions of large turbulence kinetic energy are observed to coincide with 
regions of strong vorticity, and (2) one can show that the modification of the Smagonnsky 
model which is based on the vorticity rather than the large scale strain rate, as proposed 
by Mansour et a I. to deal with regions of turbulent/ non-turbulent flow, can be shown to 
result from a production equals dissipation argument with the production term modeled on 
the rotation tensor. 

As mentioned the Smagorinsky model is based on the stress being proportional to the rate 
of strain as in a laminar flow. Lund et a 1 (1991) points out that for turbulence transport this 
relation may need to be augmented. He suggests that rotation should be included in a model 
for turbulence transport. He then proposes a subgrid-scale stress model which is based not 
only on the strain rate but on the rotation rate and products of the strain and rotation rate 
tensors. Recently, Taulbee (1992) developed an improved algebraic Reynolds stress model. 
His expression for the Reynolds stress is linear in the rate of strain but contains higher-order 
terms involving combinations of the rotation and strain tensors. Taulbee showed that the 
improved algebraic model gives better predictions that the previous algebraic Reynolds stress 
models for the simple flows he considered. 

The turbulence Prandtl number, < 7 *,, is taken as unity and Co is a model constant chosen as 
0.5. We shall now consider closure of the thermochemical equations. 


3 Thermochemical Formulation 


We consider a two-species reaction of the type A -f B — ► Products in an isothermal, incom- 
pressible flow. The species A mass fraction, Y A , conservation equation, with D constant, 


is. 


dY A , du t Y A „ , • 

+ -x = P o 


dt ' dxi dxidxi 

and a similar equation for species B. Filtering the above equation, 


( 17 ) 


Let, 


dY A du t Y A 8 2 Y A 

dt dxi dx^Xi A 

Y a = Ya + Y' a , u, = TT, + u' 


(18) 


(19) 
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then, 


( 20 ) 


ot a m Y± + g IS ± + ^aL t s1Ya r~ 

di dx{ dxidxi 


where u~ Y A is computable from the resolved field and, 


it, 


1Y' a + u\Y a + uJA' = X>T 




(21) 


So then, 


dv A , ^ (,„ , ^ 


( 22 ) 


Here Pj = vt/Sct where we choose Scj = 1.0. 


4 Chemical Source Formulation 


The chemical source term appearing in the species continuity equation is of the form 


U A = -k;Y A Y B 


( 23 ) 


where kj is the rate constant and is quantified by the Damkohler number, Da , is defined as: 


Da 


kf Y Aoo 

A U/X 


( 21 ) 


where Y Aoo is the freestream species concentration, Af7 is the freestream velocity difference 
and A is the vorticity thickness. The filtered non-dimensionalized source term is, 


ZZ=-Da(Y A Y B + Ym) ( 25 ) 

An alternative approach is to replace the spatial average with a probability average, that is, 

ZT A = -Daj u A V AB W,r)dtP'dr (26) 

where ip 1 and ip" are the sample space variables for species A and B. Thus if we knew the form 
of the joint PDF of species A and B we would be able to determine the space filtered chemical 
source term exactly. Based on our experiences as indicated in the previous section, here we 
employ a member of the Pearson family of distributions, namely the Dirichlet distribution, 
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for the joint PDF of A and B. We note that this assumption implies a univariate Beta density 
for the marginal species distributions. 

The Dirichlet or joint Beta distribution for the two species A and B is defined as follows, 


VabWW) 


r(pi + P2 + P3) 

r(pi)r(p 2 )r(p 3 ) 


(t/> , ) p, - 1 (j/>") P2-1 (i - 4 >' - ") P3_1 


(27) 


where, 


4>' > 0 , Ip" > 0 , Ip' + Ip" < 1 , Pl,P2,P3 > 0 


(28) 


The parameters pi , p 2 , p 3 are determined from the knowledge of any three of the following five 
quantities: Y\, Yb, subgrid scale species A variance, Ff , subgrid scale species B variance, 
Y§, or the covariance of the two species, Y^Yg. The chemical source term should depend on 
all of the above turbochemical quantities. We shall discuss this later. If we select the two 
species mean values along with the species covariance we can determine the three parameters. 
They are coupled through: 

Pi = -Y a S (29) 

p 2 = -Y b S (30) 

P3 = {y~A + n-\)S ( 31 ) 


where, 

s = (32) 


With this assumed form of the joint PDF, one can simply integrate the chemical source term 
in order to obtain its large scale component. This results in closed form expression for the 
source term, 


u>a = —Da 


P 1 P 2 

(u T 1 )<x 


(33) 


where a = pi -fi p 2 + p 3 . Substituting the definitions from Equations (26-29) into Equation 
(30) it can be seen that the source term is consistent with Equation (22). This serves to 
highlight the difference between mean chemistry and accounting for the subgrid fluctuations. 
The values of Ya and Yb are computable from their respective transport equations. In order 
to ascertain the subgrid species covariance Y^Yg we must solve an additional transport 
equation. This equation can be derived in a straighforward manner and its modeling is well 
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known in turbulence literature (Pope, 1979): 


pym 

Dt 


= _ (z> + p,) 


dYM 

dxi 


4 - C\Dt 


' dY A dY B 

^ dxi dx{ 


ki 


-CtYW-^ + UAYb+UBYl (34) 


where C\ and C 2 are adjustable constants. 

The chemical source terms in this equation can be likewise modeled using the assumed PDF. 
After integration their form becomes, 


CjbY'a = -Da 


(p 1 + ])plP 2 y- P 1 P 2 \ 

(a + 2)(a + l)a A (a+\)a) 


(35) 


U) A ^ g — Dd 


(P 2 + l)piP 2 _ y- P 1 P 2 \ 
(a + 2)(a + l)a e (a+l)a/ 


(36) 


The TKE and the subgrid scalar covariance are initialized with their respective time zero 
filtered DNS quantities for comparison purposes. Thus the hybrid Smagorinsky/PDF subgrid 
scale closure is complete. 


5 Numerical Formulation and Problem Description 

In order to solve the set of Equations (10,19,32), boundary and initial conditions are needed. 
In the homogeneous flow case the 2D simulations are performed with a pseudospectral col- 
location alogorithm utilizing 256 2 collocation points. The simulations are similar to those 
described in Madnia et a/. (1992) and the reader is refered to this for more details. In 
the shear flow configuration free-slip boundary conditions are employed for the cross-stream 
direction. A zero-gradient outflow boundary condition on the velocity and species fields is 
employed. The concentration fields are initialized with an error function distribution. The 
initial conditions for the mean stream wise velocity are by a hyperbolic tangent velocity pro- 
file, and the mean cross-stream velocity component is set equal to zero. In order to develop 
the vorticity rollup and pairing, which are known features of such flows, a small perturba- 
tion in the form of the most unstable mode and its subharmonics, calculated from linear 
stability analysis, is added to this mean profile. The presence of the fundamental mode in 
the shear layer produces a single vortex rollup, while the superimposed subharmonic pertur- 
bation is responsible for the second rollup in the form of the merging of two vortices. These 
phenomena have significant effects on the chemical reactions that occur within the layer. 
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The assumption of even periodicity of the flow in the y direction allows us to use pseu- 
dospectral Fourier sine/cosine expansions. However, in the streamwise direction an overall 
second-order finite difference scheme is used with a quadratic upwind differencing for the 
convection term. This eliminates the spurious oscillations associated with the dispersive 
error of the second-order scheme. The pressure is decoupled from the momentum equation 
by solving the appropriate Poisson equation. Time advancement is accomplished with the 
second-order Adams-Bashforth technique. 

The computational domain was selected to be a region bounded by (0 < x < 326) and 
(—86 < y < 86) where 6 is the vorticity thickness. For the LES runs there are 64 equally 
spaced finite difference grid points in the streamwise direction and 32 Fourier modes in the 
transverse direction. For DNS comparison runs the resolution is improved to 512 x 256 and 
for a priori analysis subsequently filtered to 64 x 32. 

A number of non-dimensional parameters characterize the flow. They are the Reynolds 
number, Re = AU8/v, based on the initial shear layer thickness, the mean velocity difference 
across the layer and the kinematic viscosity, the velocity ratio, R = AU/(U\ + f/ 2 ), and 
the Peclet number, Pe = AU8/V. Moderate values of these parameters must be chosen 
due to the computational limitations of the current simulations. It is hoped that with the 
development of the SGS model herein and future advancements these artificial limits can 
be surpassed. For this preliminary study we choose Re = Pe — 1000 and R — 0.5. The 
resolution provided here is sufficient for these simulation parameters. 


6 A priori Analysis 

The Dirichlet distribution is used to model the statistics of the species field within the 
subgrids in LES of homogeneous turbulence. Since the behavior within the subgrid can be 
assumed isotropic with a good degree of accuracy (Ferziger, 1981), the use of these models 
with their demonstrated capabilities in homogeneous flows is well justified. For the Dirichlet 
model in equilibrium chemistry, the first two moments of an appropriately defined Shvab- 
Zeldovich variable can be obtained by a properly devised SGS closure (Antonopoulos-Domis, 
1981). With the adoption of a Dirichlet density, all the higher order statistics of the reacting 
field can then be determined. In these cases, the mixture within the subgrid is usually non- 
stoichiometric, even if the initial sample within the subgrid is supplied with equivalence ratio 
of unity. In finite rate chemistry the Dirichlet distribution will be assumed for the joint PDF 
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of the two reacting species. 

In order to examine the applicability of the Dirichlet density for subgrid modeling, we have 
performed a priori evaluation of the model. These evaluations are similar to those of previous 
assessments by Erlebacher et a 1. (1987). With the construction of the DNS database for a 
2D homogeneous flow, the results are used to construct the PDF’s of the scalar variables 
within the subgrid. In Figure 1 the contours of the subgrid product mass fraction predicted 
by the Dirichlet density model are compared with those generated directly via DNS. These 
results are taken from an arbitrary plane within the three-dimensional box. The comparisons 
shown in these figures reveal a good agreement between the model predictions and the DNS 
data. This agreement provides a justification for recommending the Dirichlet density as 
a reasonable closure to account for the statistical variations within the subgrid. Due to 
a rather small sample size within the subgrid, the comparison of higher order statistical 
quantities generated by the model with those of DNS suffers from statistical errors. Also, 
the implementation of the model requires accurate input of the first two moments of the 
Shvab-Zeldovich variable. 

It has been shown previously that the Dirichlet PDF performs well in RANS modeling of 
non-homogenous flows in the limits of frozen and equilibrium chemistry (Frankel et a 1., 1992). 
With this assumed form we can obtain analytical relations for all the moments of the two 
species (Madnia et &]., 1992). We first consider the case for Da = 0. To make comparison 
with DNS results, the Dirichlet PDF is parameterized by the species mean values and the 
covariance. With this choice, in Figure 2 we display the a priori predictions and the DNS 
computations. Here we can see that the assumed PDF does very well in predicting the 
statistical behavior of the scalar values. For finite rate chemistry, consider Figure 3 for 
Da = 5. Specifically, we consider the third order moments Ya' 2 Y b and Y'^Yb ' 2 , which are 
needed in the solution to the covariance transport equation. Figure 4 shows the third order 
moments at the same Damkohler number for the case where the PDF is parameterized by 
the two means and the sum of the variances. This choice of parameters has been utilized 
by Narayan and Girimaji (1992), and requires solution of a transport equation for the sum 
of the variances, or the so called scalar energy. The differences in the statistical predictions 
between this choice and the selection of the covariance is not sufficient to warrant one choice 
over the other. Since the covariance appears naturally in the filtered source term we chose 
to use that version in this work. 

For our non-homogeneous shear flow we compare predictions obtained by employing the 
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Dirichlet distribution for equilibrium chemistry with the filtered DNS results. In Figure 5 we 
show the results for the instantaneous product thickness as a function of streamwise distance. 
The plots confirms that the Dirichlet distribution does a reasonable job of predicting the 
subgrid statistics in the limit of fast chemistry. The issue of the models validity for the case 
of finite-rate chemistry is taken up in the following section. 


7 Reacting Shear Flows 

In this section a number of simulation results are presented in order to make comparisons 
and assess the validity of the proposed SGS closure. We begin with a nonreacting flow case 
in order to ascertain whether our hydrodynamic subgrid model and our modeled species 
covariance transport equation are performing satisfactorily. This will consist of comparisons 
between two sets of simulations, LES conducted on a 64x32 grid and a DNS conducted 
on a 512x256 grid, subsequently filtered to a 64x32 grid. Next, we will consider reacting 
flows in which the performance of our hybrid Smagorinsky/PDF closure will be assessed. 
In order to see the effect of the PDF closure we shall also compare to results obtained 
using the Smagorinsky closure while neglecting the effect of subgrid scale fluctuations on 
the chemistry. These will be denoted as using the mean closure. Contour plots, transverse 
profiles and integral thickness plots of various flow-field related quantities will aid us in our 
model assessments. 

In order to assess whether the one-equation Smagorinsky model is performing adequately 
we present LES results for a nonreacting case. The filter width was chosen so as to contain 
the same physical space as the DNS results in order to make an accurate comparison. In 
Figure 6 we show transverse profiles of the streamwise velocity component at five downstream 
locations at a particular time during the simulation. The large scale velocity field is predicted 
well in the LES results. For flow visualization purposes contour plots of species A are shown 
in Figure 7. The effect of the subgrid scale model is to smear out the small scale features of 
the species field within the large scale coherent structures as evidenced in the plot. 

In Figure 8 we show contour plots of the subgrid scale TKE obtained from the (a) solution of 
the modeled transport equation and the (b) filtered DNS results. The reasonable agreement 
can be more clearly seen in Figure 9, which shows transverse profiles at select downstream 
locations. Good agreement is obtained at four out of the five downstream locations. At the 
last downstream location the TKE equation severly underpredicts the turbulence level. This 
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may be due to deficiencies in modeling the Reynolds stress in terms of the rotation tensor 
alone. Inclusion of other terms such as the strain rate and/or products of the strain and 
rotation tensors may be necessary to improve the TKE predictions. This will be discussed 
further in the next section. 

Figure 10 shows contour plots of the subgrid scale species covariance obtained from the 
modeled transport equation and the filtered DNS for the non-reacting flow case. This plot 
serves to highlight the deficiences with the modeling of the scalar flux via gradient transport. 

For the reacting flow cases considered here the Damkohler number is selected as Da = 10. 
This value is fixed for all the simulations. In Figure 11, contour plots of the product species 
from the LES and DNS are presented. From these figures we can see the LES and DNS 
results show similar large scale structures. Figure 12 shows product thickness distribution 
from the filtered DNS and the LES results obtained with and without the PDF closure. Here 
we can see that inclusion of the subgrid fluctuations gives results that predict the correct 
trend are closer to those of the filtered DNS. We note that the mean chemistry results show 
a gross overprediction of the extent of reaction when compared to the filtered DNS results. 
This highlights the importance of accounting for the subgrid species fluctuations in such 
flows. 

Figure 13 depicts contour plots of the results for the subgrid species covariance obtained 
from the filtered DNS as well as the modeled transport equation. This demonstrates the 
problem of accurately predicting the covariance. Figure 14 shows this in a more quantita- 
tive manner by displaying the average unmixedness versus streamwise distance. This poor 
agreement is in part due to the third order moments which are closed using the assumed 
Dirichlet distribution. As demonstrated in the a priori analysis section this distribution has 
difficulties with predicting the higher order moments well. The discrepancies in the covari- 
ance solution are also observed in the non-reacting case as well, indicating that there may be 
other problems with this modeled equation. These solutions were not able to be improved 
by simply adjusting the model constants. One of the major problems we feel is that the 
modeling of the scalar flux in terms of the large scale scalar gradient is inadequate in the 
context of an LES. That is, the scalar flux may depend on other terms such as the strain and 
rotation tensors. This can be deduced by developing an algebraic model for the scalar flux 
in much the same way as has been done for the Reynolds stress. This issue will be further 
discussed in the next section. 
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8 Conclusions and Extensions 


In this work, a priori and a posteriori of both homogeneous and non-homogeneous reacting 
flows are conducted in order to assess the validity of assumed PDF models as potential 
subgrid scale closures for LES of reacting flows. Specifically, a priori these studies revealed 
that if a Dirichlet distribution is assumed for the joint PDF of the two species, then for 
frozen and very fast chemistry chemistry the model performs reasonably well in predicting 
the statistical behavior of the scalar fields. For non-equilibrium chemistry, the PDF closure 
is not capable of modeling higher order moments accurately. 

LES of reacting shear flow has also been conducted employing a hybrid one-equation Smagorin- 
sky/PDF model for the statistical variations on the subgrid. Solutions for the subgrid turbu- 
lent kinetic energy and the species covariance have been obtained via solutions of modeled 
transport equations. The Dirichlet distribution has been employed to close the chemical 
source terms appearing in the thermochemical equations. The results demonstrate the effect 
of accounting for the subgrid scale species fluctuations on the product distribution within 
the shear layer. While the trends are encouraging they do highlight a number of areas where 
further work is required to improve the model. 

First, a better model for the subgrid Reynolds stress is required to obtain more accurate 
predictions for the subgrid turbulence energy levels. This might involve accounting for some 
of the additional interactions between the strain and rotation tensors as discussed by Lund 
et a 1. and Taulbee. Along the same lines, modeling of the scalar flux via gradient transport 
in the context of a time-accurate LES calculation appears unable to give accurate solu- 
tions for second-order subgrid turbulence quantities, such as the species covariance. This 
modeling deficiency will manifest itself in any potential subgrid chemistry model requiring 
such statistical information. Algebraic scalar flux models along the lines of Taulbee’s alge- 
braic Reynolds stress model would perhaps reveal some of the important terms, other than 
the scalar gradient, which would better describe the scalar-turbulence interactions. Theses 
would probably involve the strain and rotation tensors, as well as the scalar gradient, in 
some appropriate form. Finally, the discrepancies between the assumed form for the joint 
PDF of species and the actual subgrid PDF will lead to errors, as seen, when predicting the 
higher-order moments which are needed to close the equations for the second-order subgrid 
statistics. One possible approach would be instead of assuming the form of the subgrid PDF 
to solve the appropriate PDF evolution equation for the species. At this stage this lias not 
been looked at and here again some of the previous modeling assumptions may have to be 
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revised in the context of an LES. For more complex multi-species chemistry the Dirichlet 
distribution generalizes to n random variables very nicely while still providing analytic forms 
for all the joint statistics. 
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Figure Captions 


Figure 1. Contour plots of filtered product concentration. 

Figure 2. Statistical moments vs. time for Da = 0 case. 

Figure 3. Third order moments for the Da = 5 case parameterized by means and covariance. 

Figure 4. Third order moments for the Da = 5 case parameterized by means and scalar 
energy. 

Figure 5. Product thickness for equilbrium chemistry. 

Figure 6. Transverse profiles of streamwise velocity at four downstream locations (a) 55, (b) 
105, (c) 155, (d) 205. 


Figure 7. Contour plots of species A for Da = 0 case (a) LES, (b) DNS. 

Figure 8. Contour plots of subgrid scale turbulent kinetic energy (a) LES, (b) DNS. 

Figure 9. Transverse profiles of subgrid scale turbulent kinetic energy at four downstream 
locations (a) 55, (b) 105, (c) 155, (d) 205, (e) 255. 

Figure 10. Contour plots of subgrid scale species covariance for non-reacting case (a) LES, 

(b) DNS. 

Figure 11. Contour plots of product species for Da = 10 case (a) LES, (b) DNS. 

Figure 12. Product thickness for Da = 10 case. 

Figure 13. Contour plots of subgrid scale species covariance for Da = 10 case (a) LES, (b) 

DNS. 

Figure 14. Covariance thickness for Da = 10 case. 
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5 Appendix III 


This appendix provides a report, albeit incomplete, in validating the use of PF generated 
family of PDF’s in Reynolds averaging procedures. We realize that this report is incomplete 
and we apologize for that. However, we feel it is necessary to include this summary to show 
our current status in this particular aspect of our work. We are presently analyzing our 
results, and we hope to have this report completed in the near future. 

This appendix has been prepared by Mr. George Sabini, an undergraduate Research Aid, 
who is a new-comer to our group and is assisting us in this project. 



Modeling of Turbo-Chemical Fluctuations 
in a Reacting Scalar Mixing Layer 

by 

G. J. Sabini, S. H. Frankel, C. K. Madnia and P. Givi 


Abstract 


Single-point Probability Density Function (PDF) methods have proven very useful for 
modeling of turbo-fluctuations in statistical descriptions of turbulent reacting flows. In 
this work, results are presented of statistical predictions of an incompressible, turbulent, 
parallel mixing layer under the influence of a non-premixed, isothermal chemical reaction 
of the type A + B — ► Products. The self-similar region of the reacting layer is consid- 
ered. This flow configuration has been the subject of numerous previous investigations, 
and abundant experimental data are available providing the basis for appraising the per- 
formance of turbulence closures. The effects of hydrodynamics on scalar transport are 
modeled by means of conventional turbulent closures, and the influences of scalar-scalar 
fluctuations are taken into account by assumed PDF methods. Based on earlier find- 
ings, several members of the Pearson Family of PDF’s are considered: namely, a Dirichlct 
density for non-equilibrium chemistry flow, and a Beta density of the first kind for the 
flow under chemical equilibrium. The predicted results are compared with those obtained 
based on a joint Gaussian PDF, which has been employed in most previous analytical 
investigations. All the analytic results are also compared with the experimental data of 
Saetran et al (1989) and of a reacting flow and Bilger et al (1991) of a reacting flow under 
similar fluid mechanical-chemical conditions. The outcome of this comparative study in- 
dicates that the Pearson family of PDF’s are indeed superior for predicting the influence 
of turbulence on the reactant conversion rate. In particular, it is shown that the Dirichlet 
distribution if parameterized with the scalar energy, provides the most reasonable means of 
modeling the joint PDF of the scalar. The extent of agreement improves as the magnitude 
of the Damkohler number is increased. This is demonstrated by detailed comparisons of 
predicted results with laboratory data. 

Key Words: 
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Scalar Mixing Layer, Turbulent Reacting Flows, Dirichlet Distribution, Probability Density 
Function. 

Introduction: 


In our previous works (Madnia et al y 1991a, Frankel et a /., 1992) we have demonstrated 
that the Beta density provides a very good means of approximating the probability distri- 
bution of a conserved scalar quantity in both homogeneous and nonhomogeneous turbulent 
flows. Closed form analytic expressions have been obtained for the statistics of the passive 
scalar field, and with the assumption of fast chemistry, for the reactive scalar field as well. 
These relations are valid for any value of the stoichiometric coefficient. 

In this work we intend to investigate the validity of our model in predicting the statistical 
behavior of the passive and reactive scalar fields in a homogeneous mixing layer, often 
referred to as a thermal or scalar mixing layer. Results obtained with our model will be 
compared to the experimental data of Saetran et a J. (1991). In order to obtain the passive 
scalar statistics we shall follow the eddy viscosity approach detailed by Libby (1975) to 
solve the passive scalar mean and variance transport equations in this simple flow. With 
the specification of the first two moments of the conserved scalar our model is capable of 
predicting the equilibrium statistics of the reacting scalar field. 

We will then extend this work to treat nonequilibrium, that is, finite rate, chemically 
reacting flowfields. This will be done by assuming a Dirichlet distribution for the joint 
PDF of the reacting scalars. For the second order reaction under consideration here, 
this PDF is also known as the joint Beta PDF. The specification of this PDF requires 
the knowledge of first moments of both of the reactants and turbulent scalar energy Q 
(which is the sum of the scalar variances), or the reactant means and covariance. The 
information with regard to these statistics is determined from the solution of modeled 
transport equations. One of the pleasing features of the choice of the Dirichlet PDF is 
that it results in simple closed form expressions for both the chemical source term and the 
chemical source/sink terms appearing in the mean and scalar energy transport equation. 
As before all results of the investigation obtained with the PDF model will be compared 
to the experimental data of Saetran et a 1. (1991), as well as to the solution of the same 
transport equations using a Joint Gaussian PDF, which has been employed in previous 
calculations. 
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Description: 


A scalar mixing layer is formed when traditional grid generated turbulence is augmented 
with a scalar gradient. This can be in the form of either temperature, obtained as a result of 
heating half the grid (Libby, 1975) or by including two different species (Bilger, 1991). This 
simple flow configuration, because of its convenience, is of interest to study in order to gain 
insights and test modeling assumptions for more complex flows, such as nonhomogeneous 
shear layers. Experiments have been conducted by a number of researchers in order to 
study the statistical behavior of the passive scalar (Libby, 1975; LaRue and Libby, 1981; 
LaRue et a 1., 1981, Gibson et ai., 1989). Reacting scalar measurements were obtained 
by Saetran et a 1. (1989) and recently discussed by Bilger et a], (1991). Theoretical 

and modeling contributions originated with Libby (1975), Newman et ai. (1981) and 
Elghobashi and Launder (1983), and concentrated on traditiomil moment methods for 
closure. These numerous investigations have been concerned with predicting the behavior 
of the scalar mean and intensity. Particular difficulties have been noted in predicting the 
peak intensities for the passive and reacting scalars (Saetran et a L, 1989). When dealing 
with reacting flowfields the added difficulties encountered by the nonlinear chemical source 
terms make PDF methods the method of choice in such situations (Pope, 1985,1991; Givi, 
1989). 

Libby (1975) discusses the diffusion of temperature in the downstream region of a partially 
heated grid. In order to predict the passive scalar statistics in the mixing layer he utilizes 
the eddy viscosity assumption to obtain an analytic solution for the conserved scalar mean 
and to close the transport equation for the conserved scalar variance. With the grid flow 
assumptions of a homogeneous velocity field, temporally invariant statistics, high Reynolds 
number, and boundary layer assumptions he comes up with an ordinary differential equa- 
tion for the passive scalar variance. This equation has the two usual constants in front 
of the turbulent fluctuation and dissipation terms which are normally chosen to provide 
a best match to experimental results (Spalding, 1972). This linear equation is discretized 
using central differences and the resulting linear system of equations is solved using a 
tridiagonal solver. 

With the provision of the conserved scalar mean and variance, knowledge of the conserved 
scalar PDF would allow us to predict the limiting values for the reacting scalar statistics. 
Frankel et ai. (1992) (see also Madnia et ai., 1992) used a Beta density and obtained closed 
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form, analytic expression for the maximum rate of reactant conversion in a nonhomoge- 
neous shear layer. The model predictions compared favorably with those obtained from 
DNS. The derivations of the analytic expression for the statistics of the reacting scalar are 
rather involved (see Frankel, 1992). Here, we only present the final results for the second 
moment of the reacting scalar, 
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Here, B(a,f3 ) is the Beta function, and a and /3 are related to the first two moments of 
the random field (Frankel et al. , 1991). J/^or,/?) is the Incomplete Beta Function, and f s 
denotes the stoichiometric value of the Shvab- Zeldovich variable, /. For unity normalized 
concentrations at the free streams, f s = 0.5. Thus, with these formulas for reacting scalar 
mean and variance we can compare our fast chemistry predictions with the experimental 
data of Saetran et al. (1989). This will follow in the next section. 

In order to treat the nonequilibrium chemistry we employ an assumed Dirichlet distribution 
for the joint PDF of the reacting scalars. This results in analytic closed form expressions 
the chemical source and source/sink terms in the mean and scalar energy equation. The 
form of the Dirichlet distribution is 
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where 'I'' > 0, 'I'" > 0, 'P' + 'F" < 1 and Pi,P 2 ,P 3 > 0. The parameters p lt p 2 , p 3 can 
be determined from any three of the following quantities: <j4>, <B>, <.4' 2 >, <B ' 2 > , 
and <.4 ,2 > + <B' 2 >. Here we select the scalar means and the sum of the scalar 
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variances, or the reacting scalar energy, Q, and we also perform a second closure using the 
reactant means and covariance. 

The Joint Gaussian distribution used for comparison is expressed as 
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where p is the correlation coefficient, defined as <A' B'> /<B' 2 ><A' 2 >. 


Presentation of Results; 


Computations of the conserved scalar mean and variance were performed with the con- 
stants in front of the turbulent fluctuation and dissipation terms as 0.89 and 5.72, respec- 
tively. The production constant is the same as suggested by LaRue et al. (19S1) but 
the dissipation level had to be increased from their value of 2.25. This is because the 
peak intensity in the experiment of Saetran et a 1. (1989) is about 30% less than the peak 
intensity of LaRue et al.. Figures 1 and 2 show the the conserved scalar mean and stan- 
dard deviation plotted against the similarity variable 77, normalized by the mixing layer 
thickness 8. 8 is the distance between the points where the mean is 0.1 and 0.9 (Saetran 
et al., 1989). Excellent agreement with the experimental data is observed. With the first 
two moments all higher order moments of the conserved scalar are available from simple 
analytic expressions (Statistics, 19xx). In figures 3 and 4 the skewness and kurtosis versus 
transverse distance are shown for both the model and the experimental data. While the 
general trends are encouraging, note that local minima and maxima are not captured. 

The participating chemical species used to obtain the reacting scalar data are nitric oxide 
and ozone. Two cases were examined by Saetran et al ., a low and a high Damkohler number 
case. The experimental results for the high Damkohler number case were deemed close to 
the fast chemistry equilibrium limit and are chosen for comparison in this study. Thus, 
using the values of the conserved scalar mean and variance, as discussed, we can now assess 
our model predictions for the reacting scalar statistics. Figure 5 shows the mean reactant 
concentrations normalized by their inlet values. Also shown is the Beta density solution 
obtained from our closed form expressions for reactant conversion. The agreement is quite 
good. In figure 6 we show the standard deviations of the reactant concentrations versus 
transverse distance. Here again we note the excellent agreement with experimental data 
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in both the peak intensity level and the overall distribution. Another important quantity 
in modeling of turbulent reacting flows is the unmixedness or normalized concentration 
covariance. With the Beta density a closed form expression for the unmixedness has been 
obtained (Madnia et a 1., 1991; Frankel et a 1., 1992). Figure 7 shows a comparison of the 
unmixedness between our predictions and the experimental data with encouraging result. 

Obviously the reason for the good agreement between our model predictions and the ex- 
perimental data is the necessary agreement between our assumed Beta PDF and the ex- 
perimentally measured PDFs. Bilger et a i. (1991) provides further experimental data, 
specifically the conserved and reacting scalar PDFs. These are presented in figure 8. The 
agreement between the measured PDFs and the Beta PDF provides the rationale for our 
previously discussed statistical concordance. Due to the excellent agreement and attractive 
simplicity of our closed form expressions for the entire statistical behavior of the scalar 
mixing layer field, and previous agreements in both homogeneous and non-homogeneous 
environments, we heartily recommend usage of these formulas for predictions of the limiting 
bounds of the reactant conversion under warranted conditions. 


Of greater interest is the prediction of species concentration for finite rates of reactant 
conversion. With the assumed Dirichlet distribution, analytic expressions for the statistics 
of the scalar field allow closure of the chemical source terms in the mean and scalar energy 
transport equations. The three parameters of the Dirichlet distribution are p\ =<<4>a, 
P2 =< B > a, and P 3 = (1— <A> — <B>)a, where a depends on whether the third 
specified quantity is the species covariance or the scalar energy. With the scalar energy as 
the specified quantity, 

- <^> (1- <A>)+ <B> (1- <B>) _ 

° <A' 2 > + <B' 2 > 

and specifying the covariance, 
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In a scalar mixing layer with homogenous turbulence and the assumptions of stationary 
flow, self similarity, and high Reynolds number, the transport equation of Q is 
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where 77 is the similarity transformation variable, X is the downstream location normalized 
by the turbulence grid mesh size, and a and 7 are the constants in front of the turbulent 
fluctuation and dissipation terms. The analogous equation for the species covariance is 


d 2 <c4'ff'> 7 7 d<AW> 0 d<A> d<B> 

a dr] 2 +2 dr] + “ dr] dr] 


— 7 <A' B'> + X(u)a.A' + ugB 1 ) = 0 


The derivatives were approximated with second order central finite differences at fifty 
gridpoints between the boundaries of rj = 4 and 77 = —4. The statistical closures in terms 
of the scalar means and energy result in a system of coupled non-linear equations which 
were solved with a Newton- Rapson scheme with a convergence criterion of 1 x 10 -5 . 

For purposes of comparison, in addition to the Dirichlet and the Joint Gaussian models, 
also used is a mean-chemistry model, in which fluctuations are ignored in the chemical 
source terms of the transport equations The solutions of all the closure models for the 
mean reactant concentrations versus the normalized similarity variable, are shown in 
figures 9 and 10, for the the low and high Damkohler number cases, respectively. In order 
to make a better evaluation of the merits of each model than is possible with the mean 
concentration profiles, presented in figures 11 and 12 (for the same respective case) are 
the mean product profiles. The Dirichlet model parameterized with scalar energy provides 
the most accurate measure of peak product generation, while the mean chemistry model 
overpredicts the reactant conversion (especially in the center of the mixing layer, where 
the reaction occurs most vigorously). 

Figures 13 and 14 show the total scalar energy of the high and low Damkohler number 
cases for all models. The Dirichlet model parameterized with scalar energy captures the 
reduction of scalar energy in the center of the mixing layer (where the bulk of the reaction 
occurs), whereas the other models do not. The negative of the unmixedness (in this case 
equal to the covariance, because the intial reactant concentrations are normalized to unity) 
are shown for Da=0.3 and Da=1.81, respectively, in figures 15 and 16. 

Bilger (1991) provides experimental data of the skewness and kurtosis of mixture fraction 
for a variety of Damkohler numbers. For all the cases presented, it appears that these 
moments are invariant with Damkohler number. Figures 17 and 18 present the skewness, 
and 18 and 19 the kurtosis, of the two cases of finite rate chemistry (again, Da-0.3 and 


7 



Da=1.81). In both mixture fraction moments, the higher Damkbhler number case better 
reproduces the experimental data. 

Again the reason for the favorable results of our model is the agreement between our 
assumed PDFs and the measured PDFs. Bilger provides experimental data of the joint 
PDFs of the concentrations of the two reactants. Joint PDFs of the species concentrations 
according to the Dirichlet model are presented in figures 21 and 22. 

Figures 23 and 24 show the contribution to the scalar energy transport equation of the 
dissipation, production and chemical terms for the two cases of Damkohler number. In a 
recent paper, Baurle et a I. (1992) use an assumed Dirichlet distribution for concentration in 
supersonic turbulent combustion. They find that the chemical source acts as a sink which 
effectively reduces the total scalar energy, and is not consistent with the experimental data. 
We can draw the same conclusion from our results. For the high Damkbhler number, 
the chemical term is a greater sink than the low Damkbhler case, almost matching the 
dissipation term at its peak. Hence the better match in peak intensities of total scalar 
energy for the low Damkbhler number case. 

As a result of the encouraging agreement for predictions of finite rates of reactant con- 
version in the scalar mixing layer and of mixture fraction moments with using the scalar 
energy-parameterized Dirichlet distribution, we extend our recommendation, in the ab- 
sence of better alternatives, to this model. 
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